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Summary 


in this grant period, the focus has been on the effects 
of thermo-chemical nonequilibrium in low-density gases, 
and on interactions between such gases and solid surfaces . 
Such conditions apply to hypersonic flows of re-entry ve- 
hicles, and to the expansion plumes of small rockets. Due 

to the nonequilibrium nature of these flows, a particle ap 
proach has been adopted. The method continues to undergo 
refinement and application to typical flows of interest . 

A number of studies have been performed for flows in 
thermo-chemioal nonequilibrium. The effects of vibra- 
tional nonequilibrium on the rate of dissociation were stud- 
ied for diatomic nitrogen. It was found that a new model 

reproduced the nonequilibrium behavior observed experi- 
mentally. This work is described in Appendix A: a reprint 
from the Physics of Fluids. A further study made compari- 
son between the chemistry models employed in particle and 
more traditional continuum computational methods. It was 
found that the 2 approaches did give good agreement under 
conditions favorable to the continuum approach: this in- 
vestigation is included in Appendix B which is a preprint 
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for an article to be published in the Journal of Thermo 
physics and Heat Transfer. Finally, a study has been com- 
pleted which considers the effects of using particle meth- 
ods to predict the radiative emission produced in strong 
shock-waves in air. In comparisons of particle and contin- 
uum results, the total radiative intensity was found to be 
the same, although details in the profiles showed signif- 
icant differences . This work was reported at the AIAA 23rd 
Plasmadynamics Conference in Nashville, July 1992 and is 
included as Appendix C. 

The particle method was used for analysis of aerody 
namic torques for stability of the Magellan spacecraft en- 
tering the atmosphere of Venus at non-zero pitch and yaw at 
hypersonic, rarefied conditions. This work is submitted 
for presentation at the AIAA Atmospheric Flight Mechanics 
Conference, August 1993 (see Appendix D) . New models were 
developed for computing vehicle surface temperatures di 
rectly in the particle simulation which account for sur- 
face radiation, heat capacity and conductivity of the ma- 
terial, and the transient nature of the aerodynamic heat- 
ing pulse during aerobraking. This work was presented at 
18th International Symposium on Rarefied Gas Dynamics 
in Vancouver, Canada, July 1992 (see Appendix E) . A para- 
metric study was initiated which determines the accuracy 
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of particle flow solutions over a wide range of Mach num- 
ber, Knudsen number, and wall temperature for Hard-sphere 
and Maxwell molecules, depending upon the grid resolution 
and flow domain employed in the simulation. This work is 
submitted for presentation at the AIAA 28th Thermophysics 
Conference, July 1993 (see Appendix F) . 
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APPENDIX A 


Analysis of vibration-dissociation-recombination processes behind strong 
shock waves of nitrogen 

lain D. Boyd a) 

Eloret Institute, 3788 Fabian Way, Palo Alto. California 94303 
(Received 9 April 1991; accepted 23 July 1991) 

Computations are presented for the relaxation zone behind strong, one-dimensional shock 
waves of nitrogen. The numerical results are compared with existing experimental data. It is 
indicated that the derivation of chemical rate coefficients must account for the degree of 
vibrational nonequilibrium in the flow. A nonequilibrium chemistry model is employed 
together with equilibrium rate data to compute successfully the flow in several different 
nitrogen shock waves. The analysis is performed with the direct simulation Monte Carlo 
method (DSMC). The DSMC code is vectorized for efficient use on a supercomputer. The 
code simulates translational, rotational, and vibrational energy exchange, and dissociative and 
recombinative chemical reactions. A new model is proposed for the treatment of three-bodv 
recombinative collisions in the DSMC technique, which usually simulates binary collision 
events. The new formulation represents improvement over previous models in that it can be 
employed with a wide range of chemical rate data, does not introduce into the flow field 
troublesome pairs of atoms that may recombine upon further collision (pseudoparticles), and 
is compatible with the vectorized code. 


1. INTRODUCTION 

The direct simulation Monte Carlo method (DSMC) 
developed extensively by Bird ( Ref. 1 , Chap. 7 ) is a powerful 
technique for simulating flows of rarefied gases. The method 
simulates the gas flow as a collection of model particles that 
move through physical space undergoing the intermolecular 
collisions and boundary interactions appropriate to the local 
flow conditions. Collisions are simulated on a statistical 
rather than deterministic basis, and local conditions are cap- 
tured by dividing physical space into a network ot cells. The 
numerical expense of the technique is directly proportional 
to the density of the flow that has restricted its application to 
flows in the transitional regime lying between continuum 
and free molecular. 

Because of the relative computational simplicity ot the 
one-dimensional normal shock wave, this problem provides 
an opportunity to extend the use of the DSMC method into 
the continuum regime. Previously, DSMC has been em- 
ployed in normal shock waves to study translational, rota- 
tional, 3 and vibrational nonequilibrium. 4 Earlier, the com- 
putation of dissociating nitrogen shock waves was 
performed by Bird. 5 However, this study involved a finite- 
difference type approach to the problem in which constant 
volume relaxation was performed in each computational 
cell. By comparing the simulated results with the corre- 
sponding finite difference expressions, the values of tem- 
perature, density, and velocity were adjusted in each cell 
over each time step. This rather cumbersome procedure was 
adopted because of the severe computational penalty asso- 
ciated with performing a full DSMC calculation of the flow. 
One of the purposes of this paper is to report upon the devel- 
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opment of a DSMC code that simulates chemical reactions 
and is highly vectorized for efficient performance on super- 
computers. The speed-up attained through vectorization al- 
lows the proper computation of normal shock waves with 
the DSMC technique. For compatibility with the vectorized 
procedures, a new model for simulating three-body recombi- 
nation reactions in the DSMC technique is introduced. The 
new model is described, and will be shown to have several 
advantages over the previous recombination procedures. 
The flow conditions of the shock waves computed in this 
study have been chosen to correspond to some of those inves- 
tigated previously in the literature. 

II. VECTORIZED IMPLEMENTATION OF PHYSICAL 
MODELS 

As described by Bird ( Ref. 1 , Chap. 7 ) , the DSMC algo- 
rithm consists of moving the particles, sorting the particles, 
colliding the particles, and sampling the properties of the 
particles. A vectorized implementation of the DSMC algo- 
rithm for a single-species gas has been described by Boyd. 6 
The performance attained with this implementation is equal 
to that first reported for such particle methods by Baganoff 
and McDonald. 7 For the computation of reacting shock 
waves, this implementation must be extended to simulate a 
mixture of gases, and to include the additional collision 
events of vibrational energy exchange, and dissociative and 
recombinative chemical reactions. 

Vectorized implementations of the tasks performed in 
the DSMC algorithm are provided by Boyd.* Of these, only 
the implementation for the sampling of the particle proper- 
ties must be revised for the simulation of a multispecies gas 
mixture. In the new implementation, a loop is made over the 
total number of simulated particles. For each of these parti- 
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cles, the current computational cell and species numbers are 
identified. Then, several counters are updated for that par- 
ticular cell and species. These counters consist of the individ- 
ual velocity components, the squares of the velocity compo- 
nents ( these two quantities are employed in the evaluation of 
translational temperature), the rotational and vibrational 
energies, and the number of particles sampled. If the parti- 
cles are sufficiently well mixed within the simulation, then 
vectorization of this implementation works perfectly. How- 
ever, should two particles of the same species and cell num- 
ber occur over the same vector length, then only the proper- 
ties of the latter molecule will be recorded. While this means 
that the information of the first particle is lost, this aspect 
does not lead to any error. The number density in each cell is 
recorded accurately as the sorting algorithm produces this 
information. 

The vectorized shock code simulates a number of differ- 
ent collisional phenomena: translational, rotational, and vi- 
brational energy exchange; dissociative and recombinative 
chemical reactions. The total number of collisions that ex- 
change translational energy is governed by the no-time- 
counter technique. 8 The number of these collisions that un- 
dergo rotational energy exchange is given by the expression 
developed by Boyd. 3 In a revised form, this probability is 
expressed as a function of the sum of the translational and 
rotational collision energies. This aspect permits instanta- 
neous evaluation 9 thus providing compatibility with the vec- 
torized implementation. 

The rate of exchange of energy involving the vibrational 
modes is simulated through application of a variable proba- 
bility that consists of two separate components. 10 The first is 
derived from the Millikan and White 11 experimental corre- 
lation for the vibrational relaxation time r, and is given by 


ly. Therefore the values of P u{ and P v2 must be obtained by 
averaging over all collisions. In the vectorized implementa- 
tion, an average vibrational exchange probability is assigned 
to each computational cell. The probability is updated by 
employing the relative velocity of every collision that occurs 
in each cell over each time step. This process is performed 
after the pairing of collision pairs, but before computation of 
the collision mechanics, and is fully vectorized. 

The mechanisms of both rotational and vibrational en- 
ergy exchange follow the usual model of Borgnakke and Lar- 
sen, 14 and are processed efficiently using precomputed look- 
up tables. It has been observed in a recent study by Lumpkin 
etaL 15 that the mechanisms of energy transfer employed in 
the DSMC technique affect the rate of energy transfer. For 
the collision number of an internal mode, the precise rela- 
tionship between DSMC and continuum definitions is given 
in Ref. 15. It is shown that the value of the collision number 
used in DSMC will be approximately half of that determined 
experimentally and employed in a continuum or kinetic 
computation. These adjustments have been included in the 
probabilities of energy transfer to the rotational and the vi- 
brational modes. 

The chemistry model employed in the study is the vibra- 
tionaily-favored dissociation (VFD) model of Haas and 
Boyd. 13 The rate of dissociation is governed by a set of chem- 
ical rate coefficients that are reproduced in the simulation 
under conditions of equilibrium. The model also simulates 
the important physical phenomenon of coupled vibration 
dissociation. This causes a reduction in the rate of dissocia- 
tion when the vibrational energy mode is out of equilibrium 
with the translational and rotational modes, which often oc- 
curs behind strong shock waves. The dissociation probabili- 
ty takes the form 



L’l “ 

vr, 


Z| 




exp 



( 1 ) 


where v is the collision rate, g is the relative collision veloc- 
ity, and co determines the nature of the interaction potential 
in the variable hard sphere (VHS) collision model. The 
constants Z x and g* are chosen to match the experimental 
correlation for the vibrational relaxation time. For nitrogen, 
this probability becomes unrealistically high at temperatures 
greater than about 10 000 K. To simulate the relaxation time 
at such elevated temperatures, an empirical correction term 
r 2 is added to the relaxation time, where 

p c2 = 1 / vr 2 = ( 1/ Z 2 )g lt0 ‘ (2) 

The value for Z 2 for nitrogen is obtained by Haas and 
Boyd. 13 The overall exchange probability is then given as 


v(Tj + r 2 ) l/P ol 4- UP V 2 

At low temperatures, the second term tends to zero, and the 
vibrational relaxation time is given by the Millikan and 
White expression. Conservely, at high temperatures, the sec- 
ond term dominates. 

As the vibrational energy exchange probability is not a 
function of the sum of the translational and vibrational ener- 
gies, the probability should not be evaluated instantaneous- 


P d = (UZ d )[(E c -E a )*'/E?]Et (4) 

where E c is the total collision energy, E a is the activation 
energy of the reaction, and E v is the vibrational energy of the 
reacting molecule. The constant 6 determines the degree of 
vibration-dissociation coupling, and the other constants, Z d , 
and \b 2 , are determined from molecular constants, the 
chemical rate data, and the value of <b. By matching dissocia- 
tion incubation distances observed experimentally to simu- 
lated data, the appropriate value for nitrogen was estab- 
lished by Haas and Boyd 13 as 6 — 3. For compatibility with 
the VFD model and the vectorized implementation, it was 
necessary to develop a new model for simulating recombina- 
tion reactions in the DSMC technique. This model is de- 
scribed in the following section. 

The vectorized implementation of these physical models 
proceeds by constructing lists of pairs of molecules that un- 
dergo different collisional events. Initially, a list of possible 
collision candidates is formed; then a list of those that do 
collide. The pairs of molecules that collide are then subdivid- 
ed into lists of those that dissociate, those that recombine, 
those that exchange vibrational energy, those that exchange 
rotational energy, and the remainder just exchange transla- 
tional energy. The list formed of particles that experience 
each of these phenomena are then processed in computation- 
ally efficient, vectorized loops. Finally, the collision me- 
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chanics of all the molecules that collide are processed in a 
vectorized loop. A symbolic representation of the vectorized 
implementation of the collision algorithm is given in Fig. 1- 
Each of the statements represents one loop; or a set of nested 
loops in the code. Of the eleven such loops, only that entitled 
“Create List Of Pairs Of Molecules” is not fully vectorized. 
The vectorized coding provides a speed-up by a factor of 10 
in comparison with the scalar implementation of the algo- 
rithm. 

III. RECOMBINATION MODEL 

In the DSMC technique, collisions are usually consid- 
ered between just two molecules. It therefore represents a 
challenge to simulate three-body recombination reactions. 
The first such model incorporated into the DSMC technique 
was proposed by Bird. 5 In this model, a lifetime was assigned 
to a binary pair, and a three-body collision treated as a 
further collision with this intermediary, pseudoparticle. 
This scheme requires a number of assumptions to be made 
concerning the relationships between the various three-body 
collision cross sections, and the model does not lend itself 
readily to inclusion in the no-time-counter collision scheme. 
Also, under exceptional conditions, the lifetimes assigned to 
the intermediary pair can be very large, such that the pseu- 
doparticle remains in the flow field for a considerable time. 
This result is physically unrealistic, and computationally 
troublesome. An improved recombination model was pro- 
posed by Haas 16 for use in the particle method of Baganoff 
and McDonald. 7 In this model, the pseudoparticles are 
created at one time step, and then are either recombined or 
are split apart into the constitutive atoms on the following 
time step. While the scheme may be successfully incorporat- 
ed into a vectorized implementation, it still has the disadvan- 
tage of producing pseudoparticles that must be processed by 
the whole algorithm (in this case, just for one time step). 
This requirement stems from the structure of the algorithm 
proposed by Baganoff and McDonald. A further disadvan- 
tage of this model is that it requires assumptions to be made 
for the form of the molecular interaction involving the pseu- 
domolecule. The objective of the present model is to remove 
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the troublesome elements of the method developed by Haas, 
while retaining the ability to vectorize the algorithm. 

In Vincenti and Kruger, 17 it is stated that, for the re- 
combination reaction between two nitrogen atoms and any 
third body M, 

N -j_ N + M-+ N 2 + M, 
the rate of formation of N, is governed by 


^ = ( * 6 „ N /. N )« w , (6) 

dt 

where « N and n u are the number densities of atomic nitro- 
gen and the third body M and k„ is the backward ( recombi- 
nation) rate coefficient. This term is related to the forward 
(dissociation) rate coefficient k f and the equilibrium con- 
stant K c in the following way: 

k b = k f /K c = (a ) /a c )T hl ' h , (7) 

where Tis the temperature, a, and b f are the rate parameters 
associated with the forward reaction, and a c and b c are those 
associated with the equilibrium constant. The parentheses in 
Eq. (6) indicate that, in the new recombination model for 
the DSMC technique, the reaction is first treated as a binary 
collision between two atoms. This allows the derivation of a 
recombination probability in the same way as for a dissociat- 
ing reaction. To reproduce the correct rate of recombina- 
tion, the probability must then be multiplied by the number 
density n V( . The task for implementing this process in a 
DSMC simulation is to relate the rate coefficient k b to a 
probabilitv of recombination P r . This probability is ex- 
pressed as a function of the collision energy E c . When this 
probability is integrated over all collision energies, and then 
multiplied by the collision rate, it must lead to the value o 
the backward rate coefficient at a particular temperature. 
Simplification of this analysis leads to the following expres- 
sion: 






Pr(Ec)f 


(- 

\kT 



( 8 ) 


where <7 is the collision cross section for atom-atom colli- 
sions, and g is the relative velocity of collision. In the new 
model, the possibility of recombination is assessed only 
when two atoms collide. A third colliding body ( either atom 
or molecule) is then chosen at random to complete the ter- 
nary collision. As the third colliding body is chosen random- 
ly from those that exist in the cell, the recombination rate is 
achieved in the simulation by multiplying the probability by 
the total number density instead of n m . In this implementa- 
tion, the total collision energy E c is comprised of two compo- 
nents. The first, £,, is the collision energy of the two collid- 
ing atoms. The equilibrium distribution of this energy for the 
VHS intermolecular potential is 



1 

T(2-w) 





(9) 


When the two atoms come together, they form a pseudopar- 
ticle and their center of mass velocity is then collided with a 
third body to determine if recombination actually occurs. 
Bird (Ref. 1. Chap. 12) shows that the distribution function 


FIG. 1. Symbolic vectorized implementation of DSMC collision algorithm. 
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of t he center of mass velocity for the two atoms is just the 

dechon t CUy d ' Stnbuti0n of the com bined molecule. No 
th ^ arS applied t0 the coI1 jsion between the 

mndom ^ h h H PSeUd °P artlcie as th e body is chosen at 
ev E fo Th U$ lhe dlstnbut ' on function of the collision ener- 
f h y r f;‘° r the P s eudoparticle and a third body is just that for 
three degrees of freedom in translational energy, and t de- 
grees of freedom in internal energy: 




i 


F(3/2 +g/2) ikT 


1/2 - 


exp - 


/I 


\kr) 


■ C/2 - 


Xexp 


, P (. 


do 


l ( E '\ in 

T(7/2 + £ /2 — co) \kT/ 
kT ) 

hon C ,° nvemence - the form of the recombina- 

non probability is chosen as 

p ,iE c ) = (WZ r )E* c . (p) 

^“ tln§ int °, Eq l (8) and it is found 

that the follow, ng result g,ves the correct temperature de- 
pendence of the backward rate coefficient: 

X = b f ~ b c - 1/2 + o (J3) 

and 


1 


kjr 

yl/2-a, 


a f /a c 


m 1 / 2 - 


x( i 

V 2(2 — co 


co)kT 0 


rn/i + c/i-,.,) 

F ( 7/2 + £ /2 — co y) ’ 
(14) 


a W n h d e a e T m3SS ° f the a tom-atom collision, 

and a 0 is a reference cross section for the atom-atom colli- 

s,°n at temperature r 0 in the VHS collision model. The value 

Sth P ° y m C Simulation must satisfy the restnction 
hat the arguments of the gamma functions must be greater 
than zeroes is satisfied by all of the combinations of the 
lues of b f and b c found for nitrogen in the literature 
In the implementation of the recombination probability 
a multiplying factor of 100 or 1000 is usually employed 
and subsequently the recombination probability is only a^ 
sessed every 1/100 or 1/1000 times that two atoms coUid 
This procedure is adopted due to the relative infrequency of 
recombination reactions in most flows of interest, and thus 

probTem t$ E ^ nUmerically efficient mann <=r to treat the 


IV. RESULTS AND DISCUSSION 

The first test of the recombination model is performed 
under equilibrium conditions. At a constant temperature 
nd density, the model must reproduce the theoretical value 
or the average recombination probability. The results of an 
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Z?*ZY1l a,c sh r; in Fi * 2 - Th ' cond,. 
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ssxr 
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kT) 

( 10 ) 

To^btain the distribution for the total collision energy 
c ' + A** it is necessary to multiply Eqs. (9) and (10) 

onh^ 6 ™ 1 ,'^ the 11 10 ' ntegra ^ e0Vere ‘ t hsr variable. The rLult 
of this analysis is 11 


Xe.\p( - 

kji N, — N) = 7. 1 4 x 10 


llj 200/70 mVmolecule/sec, 

(15a) 

- 1.5 

Xe.xp( — 113 200/7") mVmolecule/sec, 

(15b) 

^ c(N:) = L084X 10 " ex P( - H3 200/70 molecule/m 3 . 
T . .. (15c) 

The dissociation rate coefficients are due to Bvron'* and the 
equilibrium constant is taken from Vincenti and Kruger 1 ’ it 
is shown that the model described in the previous section 

dtt^;" aVera8C Pr ° bab,llty Nation pre- 

The next test is to ensure that the dissociation and re- 
combination models allow a stat.onarv gas, initially out of 
equilibrium, to relax to the expected equilibrium state Spe 
cifying the equilibrium (constant) value of dens, ty 1 nd the 
" m ' )er !' Ure r " ,h ‘ equations are e mp "oyed 

o deiarmma the equtlibnun, atom.c mass fraction q and ha 
equilibrium temperature T: ne 


a 2 /{ \ - a) = (m ni /4p) K c 


and 


T = 


J_ r 

l-a\2 ' k 


(16a) 


(16b) 


the i 7 1S the , molecu,ar ma ^- It is assumed initially that 
the translational and rotational modes are in equilibrium 

and that the gas is purely diatomic. For the values of p = 3 

Xo ^ “ found ,hat « - w - 

pm i °?° K ' ThC Udy ° f constant volume relaxation is 

aslrvm “ “"'“T' difeenl Firs,, i, is nee 

essary to ensure tha, the dissociation and recombinatton 
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models lead to the correct equilibrium values predicted in 
Eqs. (16). Second, the effect of using different chemical rate 
coefficients is considered. The dissociation rates employed 
are those in Eqs. ( 15) and the following coefficients due to 
Kewley and Hornung: 19 


k ,(N 2 -N 2 ) = 3.82X 10 “ l r 


Xexpl — 


113 200 \ 


T 


- 


m 3 /molecule/sec, 

(17a) 


k f (N 2 -N) = i.4ixio- 4 r“ 2 - 5 

X exp ^ j m 3 /molecule/ sec . 

(17b) 

The temperature dependence of the rate coefficients in Eqs. 
(15) and (17) is significantly different. In these simulations, 
a value of <fi = 0 is employed in the dissociation model. Final- 
ly, to assess the importance of three-body reactions, the sim- 
ulation is performed without recombination reactions. 

The results of each of these investigations are given in 
Fig. 3, where the variation of a with time is shown. The solid 
curve represents the relaxation behavior using Byron s rate 
data. The chemistry models drive the mass fraction of atom- 
ic nitrogen to its expected equilibrium value. The dotted 
curve represents the simulation performed with the rate co- 
efficients defined by Kewley. With these values the simula- 
tion, again, attains the equilibrium value, but the rise in ct is 
significantly slower than for the computations that em- 
ployed the coefficients of Byron. This aspect of the simula- 
tions will be given further consideration in the computations 
behind strong shock waves. The dashed curve in Fig. 3 
shows the relaxation behavior when the recombination reac- 
tions are omitted. The value of ct will continue to rise until 
the gas consists purely of the atomic species. The effect of the 
absence of the recombination reactions only becomes appar- 
ent after the creation of a significant atomic nitrogen popula- 
tion. Further results from these three simulations are given 
in Fig. 4, where the vibrational temperature is shown. With 
Byron’s rate coefficients, the expected equilibrium tempera- 
ture is attained after a small maximum in vibrational tem- 



FIG. 3. Variation in time of the atomic mass fraction of a chemically relax- 
ing gas. 



t/t_c 


FIG. 4. Variation in time of the vibrational temperature ol a chemically 
relaxing gas. 


perature is reached. With Kewley s data, the peak vibration- 
al temperature is much higher. This effect is due to the 
reduced rate of dissociation observed in Fig. 3. However, the 
vibrational temperature does finally reach the same steady- 
state value predicted in Eq. ( 1 6b ) . By failing to include the 
recombination reactions, it is observed that the temperature 
continues to decrease as dissociation continues unabated. 
The results in Figs. 3 and 4 serve to verify the chemistry 
models employed in the DSMC technique. The need for the 
inclusion of recombination in the simulation is demonstrat- 
ed. The simulations also indicate that significant differences 
in flow properties may be computed with different rate coef- 
ficients. 

Having verified the chemistry models under test condi- 
tions, they are now applied to the relaxation zone behind 
strong shock waves of nitrogen. Three different sets of condi- 
tions are considered and are listed in Table I. The subscript 1 
indicates upsteam conditions, subscript 2 indicates the point 
immediately after the shock that is in translational and rota- 
tional equilibrium, and no subscript indicates the down- 
stream conditions. The flow parameters have been chosen to 
correspond to the conditions investigated experimentally by 
Kewley and Hornung, 19 and computationally by Bird. 5 In 
Fig. 5, the density profile is shown for case 1 . The two sets of 
computations were performed first with the rate data of 
Kewley, then with the data of Byron. In each case, the vibra- 
tion-dissociation coupling parameter 6 was set to zero. This 
leads to a minimum of coupling between these phenomena. 
It is satisfying that the computations performed with Kew- 
ley ’ s data does follow his experimental data. As may have 
been anticipated from Fig. 3, the computation that employed 


TABLE I. Experimental conditions of shock waves investigated computa- 
tionally. 


Case 

u y (km/sec) 

p t (kg/m-') 

7\(K) 

p/p i 

1 

7.31 

7.48X10--' 

25 013 

14.9 

2 

5.60 

2.86 X. 10-' 

15 130 

U.4 

3 

4.80 

4.o7 x 10 - : 

11 405 

10.0 
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Distance (cm) 

FIG. 5. Comparison of density profiles for case 1 : computed with two di ffer- 
ent sets of rate coefficients and the same chemistry model. 


Byron’s data dissociates at an elevated rate in comparison 
with the experiment. The difference is primarily in the first 
few millimeters of the flow behind the shock. 

The degree of vibrational nonequilibrium in the flow 
computed with Kevvley’s rate data is observed in Fig. 6, 
where the translational and vibrational temperatures are 
shown. For these flow conditions, the probabilities of vibra- 
tional energy exchange and dissociation just behind the 
shock are of the same order of magnitude. Hence, the cou- 
pling of the vibrational and dissociative processes is signifi- 
cant, and will limit the rate of reaction. Clearly, this phe- 
nomenon is most significant in the first few millimeters of the 
flow behind the shock. In the analysis performed by Kewley 
and Homung to derive the rate coefficients, no coupling 
between vibration and dissociation was included. Hence the 
vibrational nonequilibrium behavior had to be included in 
the data by decreasing the temperature exponents and ad- 
justing the leading constant. Alternatively, for Byron’s ex- 
periments, the shock temperatures were significantly lower, 
so that the probability of vibrational energy exchange was 
several orders of magnitude higher than the dissociation 
probability. Under these conditions, the vibrational mode 
will be fully equilibrated prior to the onset of dissociation. 
Therefore it can be stated that the temperature dependence 



FIG. 6. Comparison of temperatures for case 1: computed with Kewley’s 
rate data. 


of the rate coefficients of Byron represents the behavior in 
vibrational equilibrium. Application of the equilibrium dis- 
sociation model [6 = 0 in Eq. (4)] in the nonequilibrium 
flow therefore gives a rate of dissociation that is too high. 
There may, therefore, exist the possibility that the density 
profiles observed experimentally by Kewley and Homung 
can be reproduced with Byron’s equilibrium rate coefficients 
by applying an appropriate vibration-dissociation coupling 
model that will include the nonequilibrium effects immedi- 
ately behind the shock. 

One result of this approach is given in Fig. 7, where the 
density profiles behind the shock are again shown. By em- 
ploying Byron’s equilibrium rate data with the nonequilibri- 
um model for nitrogen [6 — 3 in Eq. (4) ], it is found that 
the experimental data is very nearly reproduced. This en- 
couraging result implies that rate data taken at moderate 
temperatures may be applicable at much higher tempera- 
tures, provided that an appropriate nonequilibrium model is 
employed. In Fig. 8, a comparison is shown of the vibration- 
al temperatures computed using Kewley’s data with 6 = 0 , 
and Byron’s data with 6 = 3 . The computations with these 
rate coefficients were also performed for the experimental 
cases 2 and 3. The comparisons of density profile are shown 
in Figs. 9 and 10. Once again, it is found that the experimen- 
tal data is reproduced with the combination of Byron’s rate 
coefficients and the nonequilibrium vibration-dissociation 
coupled model. 

Finally, the new models are applied to the flow condi- 
tions investigated spectroscopically by Sharma. 20 The flow 
was for nitrogen at an initial velocity and pressure of 6200 
m/sec and 1 Torr at room temperature. The study provided 
rotational and vibrational temperatures at two different lo- 
cations through the shock wave. The first was in the highly 
nonequilibrium region just after the shock front, and the 
second was further downstream in the equilibrium zone far 
behind the shock. For these conditions, the DSMC code was 
configured to compute the entire shock wave flow, com- 
mencing at the initial conditions. This differs from the other 
shock computations presented above which were begun at 
the point of translational-rotational equilibrium. In Fig. 1 1, 
the experimental measurements are compared with the pro- 
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FIG. 7. Comparison of density profiles for case 1: computed with two differ- 
ent sets of rate coefficients and two different chemistry models. 
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FIG. 8. Comparison of vibrational temperatures for case I • f . . , 



files of rotational and vibrational (and translational) tem- 
perature computed with the DSMC technique. Good agree- 
ment is found between the numerical and experimental data. 
Near the shock front, the DSMC simulation reproduces 
quite well the separation in the rotational and vibrational 
temperatures. Further downstream, thermal equilibrium is 
attained with values close to the measured equilibrium con- 
dition. The small d.fference in the data in the downstream 
region of the flow may be attributed to the omission of ion- 
ized reactions, radiative cooling, and the presence of impuri- 
ties m the experiment. The differences in the data in the 
nonequihbnum region must be resolved through further nu- 
merical and experimental studies. The computations were 

wrih 0r ?T USm8 eqUlllbrium rate data of Byron, together 

model thC n ° neqUlhbrium v >bration-dissoc.ation-coupled 

rem.S? ^ S ° 1Uti ° nS preSented this Studv 

reqmred about 4 h of CPU time on a Cray YMP supercom- 

The ° Ver I0 ’vollisions -ere computed for 

the 100 000 panicles in the computational domain. Without 
vectonzed implementation of most of the DSMC algo- 
rithm, the computational cost would be increased bv one 
order of magnitude on the supercomputer. Performing the 
imputations on a modem workstation would increase the 
total computational time by a further factor of at least 5 
While it can be argued that the use of such workstations is 
more cost effective (a supercomputer costs several million 


dollars) it is proposed that a turnaround time of 200 h is 
undesirable for most researchers. 

V. CONCLUSIONS 

The relaxation zone behind several strong shock waves 
MomlT 11 ! C COmputed usln S the direct simulation 

with availahl ^ excellent a S r eement was found 

with available experimental data. The intensive computa- 

■onal cost of the simulations was reduced dramatically 

ugh execution on a supercomputer of a vectorized im- 

veemnzeri 10 " f ^ a]g ° rithm ' For compatibility with the 
ectonzed implementation, a new scheme for treating three- 

reprLen™^ ' 113 ' 100 reactlons was formulated. The model 
presents an improvement on previous recombination 
schemes, and reproduces the predicted results under test 
conditions. The most important aspect of the study relates to 
the interpretation of different sets of chemical rate coeffi 
“ For shocks that exhibit significant vibrational non- 
equihbnum, the experimental density profiles mav be com- 
puted successfully in two different ways. First, bv employing 
an equilibrium dissociation model, which has no vibration 
dissociation coupling, together with the rate data denved 
from the experiment assuming vibrational equilibration 
-second, by employing a nonequilibrium, coupled vibration- 



FIG. 9. Comparison of density 
entsets of rate coefficients and 


profiles for case 2: computed with two differ* 
two different chemistry models. 
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al-dissociation model, together with rate data taken in a sep- 
arate investigation under equilibrium conditions. Clearly, 
the latter course of action is preferred, and is an important 
finding of the present investigation. 
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Evaluation of Thermochemical Models for 
Particle and Continuum Simulations of Hypersonic Flow 

Iain D. Boydf and Tahir Gok^en} 

Eloret Institute, 3788 Fabian Way 
Palo Alto, California. 

Abstract 

Computations are presented for one- dimensional, strong shock waves that are tvpical 
of those that form in front of a reentering spacecraft. The fluid mechanics and ther- 
mochemistry are modeled using two different approaches. The first employs traditional 
continuum techniques in solving the Navier-Stokes equations. The second approach em- 
ploys a particle simulation technique, the direct simulation Monte Carlo method (DSMC). 
The thermochemical models employed in these two techniques are quite different. The 
present investigation presents an evaluation of thermochemical models for nitrogen under 
hypersonic flow conditions. Four separate cases are considered that are dominated in turn 
by vibrational relaxation, weak dissociation, strong dissociation and weak ionization. In 
near-continuum, hypersonic flow, the nonequilibrium thermochemical models employed in 
continuum and particle simulations produce nearly identical solutions. Further, the two 
approaches are evaluated successfully against available experimental data for weakly and 
strongly dissociating flows. 

Introduction 

A space- vehicle passing through the earth’s atmosphere will traverse a number of 
different flow regimes. At lower altitudes, the fluid density is sufficiently large for the flow 
to be considered in thermochemical equilibrium. However, as the vehicle ascends higher 
into the atmosphere, the molecular collision rate falls, and low-density effects become 
increasingly important. 

f Research Scientist. Mailing address: NASA Ames Research Center, MS 230-2, CA 94035. 
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Continuum methods are successfully applied to flows in which the collision rate of 
the gas is sufficient to maintain Boltzmann energy distributions for the various thermal 
modes of the gas. It is not necessary that the temperatures associated with each of the 
different modes be equal, or that chemical equilibrium prevails. Particle methods, such as 
the direct simulation Monte Carlo method (DSMC), are successfully applied to flows in 
which a reduced collision rate no longer supports equilibrium energy distributions. As the 
numerical cost of this technique is proportional to the fluid density, application has mainly 
been limited to rarefied flows. 

The computation of flow properties for the flight trajectories of many space vehicles 
require the use of both continuum and particle methods mentioned above. The interface 
between the different flow regimes is therefore of great importance. Clearly, it is desir- 
able to obtain consistent results with these numerical methods in an overlapping near- 
continuum flow regime. Although the thermochemical models employed in continuum and 
particle methods are quite different, under conditions of thermochemical equilibrium they 
are expected to provide identical solutions. The relationship between the continuum and 
particle simulation under conditions of thermochemical nonequilibrium, however, has not 
been investigated thoroughly. It is therefore the purpose of the present paper to study this 
relationship by computing typical hypersonic flows with both the continuum and particle 
simulation methods. 

Evaluation of the thermochemical models is made through the computation of four 
different cases. The flow conditions in the studies are given in Table 1 and are chosen to 
examine the effects of vibrational relaxation, dissociation, and ionization. These processes 
are considered in an accumulative sense through a gradual increase in the initial enthalpy 
of the flow. The continuum and particle approaches employed in this work are briefly 
described below. 

Continuum Approach 

In the continuum formulation, a nonequilibrium eleven species gas model for air (.V 2 , 
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’ 0 ' ^ ’ ° 2 ’ N ° + ' JV+ ’ ° + ’ C ; has been implemented. The thermal state 
o the gas is described b y three temperatures: translational, rotational and vibrational- 

eectromc. The governing Navier-Stokes equations are supplemented by the equations 

account, ng for thermochemical nonequilibrium processes. The equation set consists of 

PartlaI d,fcential eqUati ° nS: «— conservation equations for species, one 

momentum equation for quasi one-dimensiona, flow, and three energ y equations. Since the 

exper, mental data is for nitrogen gas. the subset (A 2 , A', Nf, A+, ,-) of the air model rel- 
evan, to the dissociation and ionisation of nitrogen is active. The thermochemistry model 
.3 has, rally that proposed by Park. - The relaxation time for vibrational-translational 
en ergy exchange ,s taken from Millikan and White’ with Park's modification which ac- 
counts for the limning cross-section at high temperatures. Another of Park's modifications 
concern, ng the diffusive nature of vibrational relaxation is not included, wh,ch is consistent 
"-.th the current part.cle model. For vibration-dissociation coupling, the average vibra- 
tmnal energy lost or gmned due to d.ssociat.on or recombination is taken as 30 percent 
of the dissociation energy. ' The chemical reaction rates are prescribed by Park's model 

where the basic dissociation rate is assumed to be o-ovemed bv the ^ 

° uve ‘ nea dv tne geometric average of 

translational and vibrational temperatures. 

The numerical approach to solve the governing equations is fully implicit for flu.d 
dynamics and chemistry. It uses flux vector splitting for convective fluxes, and shock 
capturing. An adapt, ve grid strategy is also implemented. For the computations in this 
paper, a quasi one-dimensional code is used and a freestream of pure nitrogen is prescribed. 
The details of the numerical method can be found in Refs. 4-6. 

Particle Approach 

The particle simulation code employed in this investigation provides modeling of the 
, rotational, vibrational, and electron kinetic energy distributions. These are 
complemented through simulation of dissociative, recombinative, ionising, and exchange 
reactions. The code is vectorised for efficient execution on a Cray-YMP. Description of 
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vectorized implementation can be found in Refs. T end 8. The bound*,- condition* 

rr 90W ^ ~ end dol 

, at . 0 f ClheS The d ° WnStream Te,OCi ‘- V " ° bta,n<!d eilher ‘he continuum ceicu- 
.on* or rent liable experiment*, date. The present simuiet.ons compute the entire 

UCtUre fr ° m UPSlrea " *° d °'“ «“»*«•■ Once the shock reaches e spec- 

'J - adjUStmCntS ~ *° ^ ^ computet, one. 

gnd to maintain a steady shock position. 

The rate of energy exchange between the translational and rotational energy modes 
—ated using a probability function evaluated using the energy of each collision. > 
ames of rotational energy exchange is performed by the Borgnakke-Larsei 

S ‘ m0d6L " ^ meCh “ ,CS exchange are 

compu e using two different schemes. The first uses the B.rgnaldte-Larsen aoproach with 

a continuous vibrational energy distribution described by a fixed number of vibrational 

reedom, The second approach, due to McDonald, allows sampling of 
post-co , ,s, on vibrational energy levels from the discrete form of the Simple Harmomc 

wZlrr^ TWS “ - « — - - - he estimated for the 

thefi T f d ' U effeCtiVe - VaneS f ” aC “ rdinS l ° th ' content of 

wh . " OTe " " the Preferred 1PPr ° aCh fr ° m a Physi “‘ ^hpotnt. The manner in 

by L \ Z ^ eXChange " P6rf0rmed “ P ar “ c ‘ e simulation is shown 

y ump m e( «1. » to a ffe ct the rate of relaxation. Therefore, all the rotational and 

h . . re aXatl ° n m ° dels emplo - ved m the Particle simulations are adjusted to match 
e continuum values by the correction developed in Ref. 13, 

*• “ “■ a— a™, .a., 

™ 77“ “ ““ - " '• - — «, imMm a. 

f P yS1Cal Phen ° men0n ° f vibration-dissociation coupling. The model contains a 

ree parameter * winch controls the degree of coupling between vibrational and dissociative 
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relaxauon processes. 'It was demonstrated in Ref. 14 that, by increasing the value of 

0. It is possible to increase the dissociation incubation time in the simulation. Also in 
Ref. 14, through comparison with experimental data, the value of i, for nitrogen was 
determined assuming Borgnaklte- Larsen mechanics for vibrational energy exchange with a 
fixed value for f.. In the DSMC code, the mode, employed for the reverse recombination 
reaction appropriate to VFD is that developed by Boyd. » Al, other chemical reactions, 

1. e., ionization and exchange reactions are simulated using the steric factor developed bv 
Bird. » The inclusion of electrons in the simulation is discussed in detail in Ref. 16. 

Chemical Rate Coefficients 

The rate coefficients employed in the reactions of interest in the present study are 
0 iven in Table 2. These are described in the usual Arrhenius form: 

k(T) = a T b exp( — E a /kT) 

where a and b are empirically determined constants, £„ is the activation energy, and T is 
the controlling temperature. Three different sets of coefficients are given corresponding to 
those used: in the continuum code; in previous DSMC investigations; and in the present 
DSMC code. The values of the activation energy used in the three sets of rate data are 

unchanged for each separate reaction. Therefore, the exponential term in the Arrhenius 
form has been omitted from Table 2. 

The rate expressions employed in the continuum code are those recommended in the 
review by Parle et * > Generally, only the forward rate constants are specified. In the 
dissociation reactions, reactions 1 and 2, the controlling temperature in the continuum 
two-temperature approach is given by T a =(TT„)i. For nitrogen dissociation, the particle 
code employs the rates of Byron '• in the Vibrationally Favored Dissociation model (VFD). 

It was shown previously by Boyd* that these rates, when used w,th the VFD model, are 
capable of reproducing vibration-dissociation coupling observed at high temperatures. 

The reverse rates for each reaction are obtained by evaluating the following temperature- 
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dependent form for the equilibrium constants proposed by Park : 1 

ln(K e (T)) = Aiz + A 2 + A 3 ln(z) -f A 4 /~ 4 - A 5 fz 2 

where the A, are constants and z=10000/T. In the continuum code, the values of .4, for 
reaction 3 are obtained from Ref. 18 while those for reaction 4 are taken from Ref. 1 . 
Unfortunately, this form for the equilibrium constant is not mathematically convenient for 
implementation in the DSMC chemistry models. However, a set of reverse reaction rates 
for use m DSMC has been determined by Bird , 19 and these have been used in a number 
of studies. These reverse reaction rates are determined by calculating an equilibrium 
constant m which the exponential terms m the electronic partition functions are evaluated 
at a temperature of 11,000 K. It is possible to compute the equilibrium constants employed 
by Bird by considering the ratio of the forward and reverse rates for each reaction. This 
has been performed for reactions 3 and 4 of Table 2 . The equilibrium constant employed 
by Bird and that used in the continuum code are shown as a function of temperature for 
reaction 3 m Fig. 1 . It should be noted that the exponential term has again been omitted 
for the sake of simplicity. For this reaction, it is found that the equilibrium constant 
employed by Bird is about 2 orders of magnitude higher than the continuum expression. 

In reaction 4, the equilibrium constant used by Bird gives values which are again generally 
higher than the continuum model. 

It is to be noted that the goal of the present study is to evaluate differences in the 
chemical models employed in the continuum and particle methods. To limit the number 
of factors involved in our comparisons, it is the aim to maintain consistency between the 
relaxation rates employed in the solution techniques. Therefore, a form for the equilibrium 
constant which takes the traditional Arrhenius form is curve-fitted as a function of tem- 
perature to Park’s expression. The limitation on the Arrhenius form which may be used 
conveniently m Bird’s expression for the probability of chemical reaction 15 is discussed 
by Boyd and Stark . 21 The curve fit for reaction 3 is also shown in Fig. 1 . The resulting 
rate constants for the reverse direction for reactions 3 and 4 are listed in Table 2 . Gener- 
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a lly> good agreement is obtained between the new DSMC expressions and Park’s results, 
particularly over the temperature range of interest, i.e. from 5,000 to 20,000 K. 

For reaction 5, the temperature dependent form proposed in Ref. 2 is not convenient 
for use in the DSMC chemistry models. In comparing the forward rate constants employed 
by Park and Bird for this reaction it is noted that Bird’s reaction rates are several orders 
of magnitude lower. Once again, a curve fit is made to Park’s expression in an Arrhenius 
form which may be employed in the particle chemistry models. The new form, which is 
given in Table 2, gives closer correspondence to Park’s results over the temperature range of 
interest. Further analyses have been performed which improve the correspondence between 
the chemical rates employed in continuum and particle simulation for all reactions in air 
involving charged species and are reported in Ref. 16. 

Presentation of Results 

Computations are performed for four different sets of flow conditions for normal shock- 
waves in pure nitrogen and these are listed in Table 1 in which subscripts 1 and 2 indicate 
upstream and downstream conditions, respectively. The upstream temperature is pre- 
scribed to be 300 K for all cases. The upstream density together with the length of the 
computational domain simulated are chosen such that the flows axe in the near-continuum 
regime. In other words, Knudsen number is small enough that the thickness of the shock 
wave is small compared to relaxation distance behind the shock. The different upstream 
flow conditions also provide increasing enthalpy: thus, the flow behind the shock is char- 
acterized in Case 1 by vibrational relaxation processes; in Case 2 by weak dissociation: 
in Case 3 by strong dissociation; and in Case 4 by weak ionization. The conditions in 
Cases 2 and 3 correspond to those investigated experimentally bv Kewley and Hornung. 22 
The results for each of these investigations axe described in the following sub-sections. The 
numerical parameters chosen for each DSMC computation followed the usual guidelines in 
setting the cell size less than the mean free-path and the time-step to be a fraction of the 
mean time between collisions. The cell size criterion is relaxed in regions far behind the 
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shock-front where flow gradients are less severe. In each case considered, the number of 

computational cells is 1,000 and the number of simulated particles is maintained at about 

100,000. 

Case 1: \ Ibrationally Relaxing Flow 

Density profiles for the first case investigated are shown in Fig. 2. Very good agreement 
is found between the continuum and particle simulation results. Two different DSMC com- 
putations are shown: the first corresponds to the use of the Borgnakke-Larsen approach 
(BL) for performing the mechanics of vibrational energy exchange with a constant number 
of vibrational degrees of freedom, Cv=l-6. This value corresponds closely to that eval- 
uated at the downstream equilibrium temperature. The second solution employed the 
discrete vibrational energy sampling approach for the Simple Harmonic Oscillator (SHOj 
of McDonald ' which automatically varies This is the first time that a comprehensive 
comparison is made between continuum and particle simulations for vibrational relaxation 
behind a strong shock. It is very encouraging to observe that, under near-continuum 
conditions, the two methods give such close agreement. 

The variation in translational and vibrational temperatures for this case are shown 
m Fig. 3. The particle solutions are obtained with McDonald’s variable ( v model. Once 
again, very good agreement is obtained between the two solution techniques. Temperature 
is generally a much more sensitive quantity to simulate than density. The close correspon- 
dence between the continuum and particle results indicates that the vibrational relaxation 
models of both approaches are very nearly equivalent. This comparison therefore lends 
strong support to the use in the particle simulation of the vibrational energy exchange 
probabilities developed for the DSMC method, 11 the correction term required to equate 
the continuum and particle relaxation rates, 13 and the mechanics of vibrational ener-v 

exchange. l - It should be noted that the degree of dissociation under these flow conditions 
is less than 1%. 

Case 2: Weakly Dissociating Flow 
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The second set of conditions considered has an increased enthalpy which gives rise to 
weak dissociation behind the shock. This case is of additional interest as it was studied 
experimentally by Kewley and Hornung 22 who employed interferograms to measure the 
variation in density behind strong shocks of nitrogen. The increase in enthalpy is revealed 
in the density profiles shown in Fig. 4 in which the normalized density rise reaches a value of 
about 10 at the downstream boundary. While both solutions give good agreement with the 
experimental data, it is clear that the particle solution provides the better correspondence. 
The DSMC profile is obtained with the variable £ v model. Comparison of the translational 
and vibrational temperatures computed through the shock are shown in Fig. 5. Again, 
a very good agreement is observed for the two sets of results. Considering the excellent 
agreement obtained in Fig. 5 between the continuum and particle methods, and also for the 
case of vibrational relaxation, it is concluded that the differences observed in Fig. 4 must be 
due to the dissociation models employed in each simulation technique. This indicates that 
the continuum two-temperature model gives a dissociation rate which is slightly slower 
than that of experiment and DSMC, In other words, for a weakly dissociating gas, the 
effect of vibrational relaxation on dissociation is overestimated in Park’s two- temperature 
model. 

The results for the mole fractions of molecular and atomic nitrogen are shown in 
Fig. 6. As expected from the previous comparisons, there is close agreement between 
the two numerical approaches with DSMC predicting slightly more dissociation than is 
obtained in the continuum solution. 

Case 3: Strongly Dissociating Flow 

The further increase in enthalpy for Case 3 gives rise to stronger dissociation effects. 
Once again, the flow conditions modeled match those considered experimentally by Kewley 
and Hornung. 22 The experimental profile of density behind the shock is compared with 
the computational results in Fig. 7. The comparison between the continuum solution and 
the experimental data is excellent. The DSMC profile is computed using the variable £ v 


9 



model and <p= 2. With this model configuration, the particle method provides excellent 
agreement with both the experiment and the continuum solution. It should be noted that 
this is the DSMC model configuration employed in the computations for Case 2. 

The translational and vibrational temperature profiles computed with the continuum 
and DSMC techniques are shown in Fig. 8. Generally, very good agreement is observed 
between the two. There is a noticeable difference in the peak vibrational temperatures 
computed by the two methods. This has quite significant implications for the estimation 
of radiative emission in such flows. The difference is attributable to dissociation-vibration 
coupling, i.e., how the vibrational energy distribution is affected by dissociation. This 
process is modeled quite differently in the continuum and particle approaches. These 
results indicate the need for experimental measurement of vibrational temperature profiles 
behind shock waves under conditions similar to those considered here. For completeness, 
the profiles of mole fractions of the neutral species are shown in Fig. 9. The stronger 
degree of dissociation for these flow conditions is very evident and, as expected, very good 
agreement is found between the two solutions. 

For this strongly dissociating case, it is found that Park’s two temperature model 
reproduces the experimental data very well. It is very encouraging that the two temper- 
ature model gives such a favorable comparison with the experimental data in stronglv 
dissociating flow as this is the regime for which the model has been developed. Indeed, 
the present comparison arguably provides the strongest evidence to date that, despite 
its weak theoretical basis, the two-temperature model does produce adequate simulation 
of strongly coupled vibration-dissociation processes. The present investigation is unique 
in that evaluation of the model is performed through direct comparison with experimen- 
tal measurements of a fundamental flow quantity. The model was previously calibrated 
against experimental data by Park 1 through comparison with radiation emission spectra, 
and by Candler 23 through comparison with shock stand-off distance. Due to the excellent 
comparisons between DSMC and experiment in Figs. 4 and 7. it is recommended that Mc- 
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Donald’s collision mechanics and the VFD model with 4>=2 be 
dissociation with the particle method. 


used for simulating nitrogen 


Case 4: Weakly Ionizing Flow 

The increase in enthalpy for Case 4 is sufficient to give rise to significant ionization 
effects behind the shock. In performing the DSMC computat.ons of the ionized flowfield, a 
steady shock solution is first obtained with the ionizing reactions omitted. After reaching 
this point, the ionized species are included and a further short transient phase in the simu- 
lation then allowed before sampling of flow properties is commenced. These procedures are 
adopted because the inclusion of electrons in the flowfield requires a reduction in compu- 
tational time-step by two orders of magnitude. To compute the entire flowfield with such 
a small time-step would require much larger computational resources. The comparison for 
density profiles computed with the numerical techniques is shown in Fig. 10. As with the 
previous cases, good agreement is obtained between the two solutions. The temperature 
profiles computed with the continuum and particle methods for the translational and vi- 
brational modes are compared in Fig. 11. The peak values for each energy mode are in 
good agreement. It is observed that the translational temperature shock computed with 
DSMC is thicker than the continuum result. This is due to the relatively low upstream 
density employed in this investigation. A more thorough analysis of such behavior will 
form the basis of future study. The computed variations in mole fractions for the neutral 
species obtained with the numerical techniques are compared in Fig, 12 and those for the 
charged species are compared in Fig. 13. The agreement which is generally obtained is 
very satisfactory. This comparison verifies that the new forms of the reverse reaction rates 
employed in the particle simulation are nearly equivalent to the expressions used in the 
continuum analysis. It should be noted that a degree of statistical scatter is exhibited by 
the DSMC results for the less abundant species. 

To assess the effect of using the new reaction rates, a particle simulation is also 
performed with the rate data used by Bird. ■» The variation in the mole fractions of the 
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charged species computed in this way are compared with the continuum results in Fig. 14. 
None of the species profiles are found to be in good agreement. With Bird’s rates, the most 
populous ion is N, + , whereas the new parttcle rate data agrees with the continuum solution 
in giving N+ as the most abundant ion. With Bird’s rate data, the mole fraction of electrons 
at the downstream boundary is about 2.5xl(r 9 whereas the continuum simulation gives a 
value of about 1.8xl<T 2 . If it is accepted that the rate coefficients provided in Refs. 1 and 2 

are the more physically realistic, these large differences observed with Bird’s older data set 
must call into question previous DSMC investigations which employed those reaction rates. 
Apparently, there has been a mistake in the reverse rate for reaction 3 used by Bird. 19 The 
temperature exponent should be listed as -0.18 instead of the value of -0.52.“ This error 
is quite serious as it was included in codes employed in a fairly large number of studies. 

Concluding Remarks 

This study was motivated by the requirement to evaluate the relationship between 
continuum and particle simulations of hypersonic flows in the near-continunm regime. The 
results obtained in the investigation have established that a close correspondence exists 
between the thermochemical nonequilibrium models employed in these solution techniques. 

In the case of vibrational nonequilibrium, the agreement between the two sets of numerical 
results validated a number of recent modeling developments for computing the rate and 
mechanics of vibrational energy exchange in the particle simulation. In the cases of weak 
and strong dissociation, both the continuum and particle models for vibration-dissociation 
coupling were successfully evaluated against experimental data. This is a most interesting 
result considering the large differences in the dissociation models employed in the two 
techniques. In the case of weakly ionized flow, it was necessary to develop new forms 
for some of the chemical rate constants for use in the particle simulation. These were 
developed so as to be nearly consistent with the continuum expressions, and also to be 
mathematically convenient for use in the particle chemistry models. The next stage in 
this continuing investigation will be evaluation of these methods for flow conditions in the 
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transition regime, i.e. at higher Knudsen numbers. In such Hows, rarefact.on affects mav 
invalidate use of the Navier-Stokes equations, and give 'rise to large differences between the 
continuum and particle simulation results. 
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Table 1. Flow conditions. 


Case 

Ui (m/s) 

Pi (kg/m 3 ) 

Pi (torr) 

L 7 2 (m/s) 

1 

4000 

1.75xl0 -3 

1.17 

541 

2 

4800 

4.67xl0 -2 

31.2 

480 

3 

7310 

7.48xl0 -3 

5.00 

496 

4 

10000 

o.OOxlO -4 

0.33 

640 
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Table 2. Leading constants in chemical rate data (m 3 /molecule/s). 



Reaction 

Continuum 2 

Particle 19 

Particle (present) 

1. 

N 2 + N 2 — + 2N + N 2 

1.16xl0~ 8 T- 1 - 6 

6.17xl0 -9 T" 1 - 6 

7.97xl0~ 13 T~°- 5 

2. 

N 2 + N -*• 2N + N 

4.98xl0 -8 T -1 - 6 

l.SoxlO -8 T -1 - 6 

7.14xl0 -8 T" 1 - 5 

3a. 

N" r + N 2 — + N 4* N2 4 ” 

1.66xl0 -18 T 0 - 5 

1.67xl0 -17 T -018 

1.66xl0 -18 T°- 5 

3b. 

N + N 2 + ->• N + + N 2 

See Ref. 1 

2.37xl0~ is T" 0 - 52 

2.34xlCT 14 T -0 - 61 

4a. 

N + N-» N 2 + + e" 

7.31xl0 -23 T 1 - 5 

2.98xlO -20 T°- 77 

7.31xl0~ 23 T 1 - 3 

4b. 

N 2 + + e" N + N 

See Ref. 1 

8.88xl0 _1 ° T- 1 - 23 

1.57xl0~ 17 T 0 ' 85 

5. 

N + e~ — + N + + 2e _ 

4.15xl0 4 T" 3 ' 82 

l.OOxlO -14 

S.SlxlO -8 T _1 -° 
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chemical rate data. 


Miszi/gt- 

PRECED'.NG PAGE -etAWt NOT FILMED 


19 




fy. I ■ ^ 







^ & a s' 

^ v 




Tt (continuum) 

Tv (continuum) 

A Tt (particle) 

V Tv (particle) 

! 1 I 


Distsnce (cm) 


. S . Sond ^ ^oh 





0.0 0.1 0.2 0.3 

Distance (cm) 







0.0 0.1 0.2 
Distance (cm) 






— N2 (continuum) 

N (continuum) 

N2 (particle) 

1 N (particle) 

00-0 OOCC.00.0-C.CC 

J 1 I 1 

.1 0.2 0.3 
Distance (cm) 


• 4 . y~- £oice^ 




Experiment 
Continuum 
Particle (SHO, 6*2} 


0.1 0.2 
Distance (cm) 


•n 

O’ 7 







F/jX 






£ qL-CjZs^ 


f: v “7. 


k 






°? 

o 


60 



Distance (cm) 


fk.il. U 





Distance (cm) 


12. £okcJ2.~^ 




Mole Fraction 



Distance (cm) 


73. 


5 . 





Distance (cm) 


Fj . ^ <£x>)4 V- £ 



APPENDIX C 




AIAA-92-2971 

Decoupled Predictions of 
Radiative Heating in Air 
Using a Particle Simulation Method 

I. D. Boyd and E. E. Whiting 
Eloret Institute 
Palo Alto, CA 94303 . 


A1AA 23rd 

Plasmadynamics & Lasers Conference 

July 6-8, 1992 / Nashville, TN 


For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics 
370 L'Enfam Promenade, S.W., Washington, D.C. 20024 


DECOUPLED PREDICTIONS OF RADIATIVE HEATING IN AIR 
USING A PARTICLE SIMULATION METHOD 

Iain D. Boyd* and Ellis E. Whiting" 

Eloret Institute, 3788 Fabian Way 
Palo Alto, California. 


Abstract 

The radiative emission along the stagnation stream- 
line and the radiative heating at the stagnation point 
of a blunt-nosed vehicle entering the Earth’s atmo- 
sphere at hypersonic speed are estimated using a par- 
ticle simulation technique with decoupled radiation. 
The fluid mechanics of the weakly ionized flow is 
computed using the direct simulation Monte Carlo 
method, DSMC. Analysis of radiation is decoupled 
from the flowfleld and is estimated using NEQAIR, 
a computer code written by Park. The results are 
compared with previous DSMC computations which 
employed a simplified, coupled radiation model. The 
effects of recent advances in modeling relaxation and 
dissociative behavior with the DSMC technique are as- 
sessed in terms of radiative emission. It is found that 
the introduction of the new models decreases the'pre- 
dicted total radiative heating at the stagnation point 
of the vehicle by a factor of 15. The DSMC approach 
is also compared with a continuum flow model: in each 
case prediction of radiation is decoupled from the flow- 
field. Similar sets of vibrational and chemical relax- 
ation rates are employed in these simulations. Despite 
large differences in the computed flowfleld, which ex- 
hibits strong thermo-chemical nonequilibrium, the to- 
tal predicted radiative heating estimates agree within 
a factor of 2. 

Introduction 

The motivation for the present study arises from 
the requirement for accurate radiation estimates for 
hypersonic flight vehicles. These are necessary for ad- 
equate thermal protection of the spacecraft during en- 
try into the Earth’s atmosphere. A previous study for 
the Aeroassist Flight Experiment (AFE) has been re- 

* Research Scientist. Mailing address: 

NASA Ames Research Center, MS 230-2, CA 94035. 

Thij paper i n declared a work of the U.S. Government and is not 
subject to copyright protection in the United States. 


ported by Whiting and Park 1 in which flowfleld data 
were obtained using the Navier-Stokes equations at the 
lower altitudes traversed by the AFE during its sweep 
through the upper atmosphere. However, at higher 
altitudes, numerical difficulties were found in obtain- 
ing solutions. These difficulties are presumed to arise 
from the failure of the Navier-Stokes equations when 
the Knudsen number of the flow is too high. 

The primary aim of the present study is to obtain 
radiation estimates using a particle simulation method 
for low-density, hypersonic flows in which strong effects 
due to thermo-chemical nonequilibrium are present 
and to compare the results with those from Ref. 1 for 
continuum flow. This is accomplished using the di- 
rect simulation Monte Carlo method (DSMC) to pre- 
dict the one-dimensional flowfleld along the stagnation 
streamline. The radiative emission is then computed 
using the computer code NEQAIR 2 with the radiation 
decoupled from the flowfleld solution. 

This is the first time that radiative heating has been 
estimated from DSMC computations using the ap- 
proach in which radiation is decoupled from the flow- 
field solution. A new particle simulation code has been 
developed for the study. The code contains many re- 
cent developments which have improved considerably 
the modeling of thermo-chemical relaxation with the 
DSMC technique. 

The new code is first evaluated through comparison 
against previously published particle computations of 
hypersonic, low-density, radiating flow. These previ- 
ous DSMC computations employed simplistic thermo- 
chemistry models, and the effect of the introduction of 
the improved models is assessed. Additionally, com- 
parison is made between the new code and continuum 
techniques for flow conditions corresponding to a tra- 
jectory point of the AFE vehicle. The continuum code 



2 AIAA 92-2971 


is SPRAP (Stagnation Point Radiation Program) writ- 
ten by Park 3 which solves the Navier-Stokes equations. 
Radiative emission for these continuum calculations is 
also decoupled from the flowfield solution. 

Description of the Particle Simulation 

The particle simulation code for computing one- 
dimensional shock waves in air is highly vectorized, 
and has been extended from the implementation for a 
reacting gas described by Boyd 4 to include electrons 
and ionizing reactions. An eleven species (N 2 , N, 0 2 , 
0, NO, N 2 +, N+, 0 2 + , 0+, N0+, E") real air model 
is employed. The code contains many recent model- 
ing developments which have improved the ability of 
the DSMC technique to simulate thermochemical re- 
laxation phenomena. 

The code simulates rotational 5 and vibrational 6 en- 
ergy exchange using probability functions which are 
evaluated using the energy of each collision. These 
represent an improvement on previous models in which 
constant probabilities were applied. The code employs 
the Vibrationally-Favored Dissociation model of Haas 
and Boyd 7 and a corresponding recombination model. 4 
From the DSMC simulation, it is assumed that the vi- 
brational and electronic temperatures of the gas are 
equal for the purposes of computing the radiative emis- 
sion. Therefore, it is expected that the improved phys- 
ical models employed in the DSMC code may affect 
significantly the total radiative heating predicted. 

The introduction of electrons into the code requires 
special consideration. Due to their very low mass, 
the collision rates associated with electrons are about 
two orders of magnitude higher than those associated 
with the heavy species which occur in air. This re- 
quires the use of a numerical time-step which is two 
orders of magnitude smaller than would be employed 
were the electrons absent from the flowfield. In the 
present implementation, this problem is solved by us- 
ing a time-step based on the heavy species to move the 
particles through the flowfield, and then subdividing 
the time-step into one hundred sub-steps to perform 
the collisions. Charge-neutrality is enforced through- 
out the flowfield by forcing each electron to remain 
in the same computational cell as the ion with which 
it was initially formed. While this is physically un- 
realistic, this approach should not affect the predic- 
tion of radiative emission to any great extent. Also, 


any charged particles which collide with the cold sur- 
face of the vehicle are assumed to recombine to the 
neutral species. The special considerations described 
here for simulating charged species add a significant 
numerical overhead to the calculations (although the 
DSMC code is still highly vectorized). This overhead 
is minimized by initially running the DSMC code with- 
out the charged species until a steady state is reached. 
The ionizition and charge-exchange reactions are then 
turned on, and the system is allowed to reach a new 
steady state, before sampling of flow variables begins. 

Two different sets of chemical rate data are em- 
ployed in the DSMC code. The first corresponds to 
that used by Moss et a/..® This set was implemented 
without coupled vibration-dissociation. In addition, 
for each reaction, the reverse rate expressions were de- 
termined by evaluating the electronic partition func- 
tions in the equilibrium constants at a fixed tempera- 
ture. The second set of rate constants is based on those 
employed by Whiting and Park in the continuum code 
SPRAP. In their analysis, the reverse reaction rates are 
obtained by evaluating temperature-dependent curve 
fits for the equilibrium constant. The curve fits take 
the following form proposed by Park 3 : 

K e (T) = exp(Ai + A 2 ln{z) + A 3 z + A 4 z 2 + A 5 z 3 ) 

where the A{ are constants and z=10000/T. Unfortu- 
nately, this form for the equilibrium constant is not 
mathematically convenient for implementation in the 
DSMC chemistry models. 

The goal of this part of the present study is to evalu- 
ate differences in the chemical models employed in the 
continuum and particle methods. To limit the num- 
ber of factors involved in these comparisons, it is the 
aim to maintain consistency between the relaxation 
rates employed in the solution techniques. Therefore, 
a form for the equilibrium constant which takes the 
traditional Arrhenius form is fit as a function of tem- 
perature to Park’s expression. This approach has been 
very successfully applied by Boyd and Gok^en 9 to the 
ionizing reactions for N 2 . By achieving good corre- 
spondence between the two sets of chemical rate con- 
stants, very good agreement was obtained for flow so- 
lutions to hypersonic shock-waves in N 2 using parti- 
cle and continuum methods. For air, good agreement 
is generally obtained between the new DSMC expres- 
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sions and Park’s results over the temperature range of 
interest, i.e. from 10,000 to 20,000 K. 

Description of NEQAIR 

The code employed to estimate the radiative emis- 
sion from the flowfield solution is called NEQAIR and 
was developed originally by Park. 2 NEQAIR calcu- 
lates the nonequilibrium electronic excitation of the 
species in the flow and the radiation emitted by those 
species. A principal assumption made in the code is 
that the quasi-steady state (QSS) model may be used 
to determine the population of each electronic state of 
a species after the species density is found from the 
reacting-flow solution. 

The QSS model assumes that the individual rates 
of populating and depopulating each electronic state 
are much faster than the rate of change of the pop- 
ulation of the state itself. Thus, the population of a 
state is given by the difference of two large numbers 
that are about equal. This assumption is not valid in 
the flow immediately behind the shock wave where the 
populations of the excited states are low and increas- 
ing rapidly. However, as the excited states become 
populated, the species begin to emit radiation, and by 
the time the radiation reaches a significant level, the 
assumption is reasonably valid. 

As the QSS model is applied after the chemistry 
portion of the calculation is completed, the electronic 
excitation processes are assumed to be independent of 
the chemical reactions. This simplifies the calculation 
enormously. The QSS model allows a separate effective 
electronic temperature to be defined for each excited 
electronic state. These temperatures are defined by 
comparing the population of an excited state to the 
population of the ground state. At each flowfield data 
point, the QSS model implemented in NEQAIR em- 
ploys the translational and vibrational temperatures 
from the DSMC or continuum simulation. In addi- 
tion, the QSS model assigns many temperatures for 
electronic excitation (one for each excited electronic 
state of every species). 

In the estimation of radiative emission from the 
flowfield solution, NEQAIR is executed over 12 differ- 
ent spectral regions which are listed in Table 2. These 
have been determined by Whiting and Park 1 to be 
those which make significant contributions to the total 
radiative heating for hypersonic flows of air. The fol- 


lowing molecular bands have been included: N 3 + (1-), 
N 2 (1+), N 2 (2+), N 2 (BH1), NO/?, N0 7 , 0 2 (SR). ’ 

The ultimate test of such codes is how well the re- 
sults compare with experiment. Park 10 has compared 
calculated radiative results with shock-tube, ballistic- 
range, and Earth-entry data covering a wide range of 
flight conditions and finds good agreement, generally 
within a factor of two. 

Results and Discussion 

This section is divided into three sub-sections. In 
the first, comparison is made between the new code 
and DSMC computations reported previously by Moss 
ti a/.. 8 A number of different comparisons are made 
to examine the ability of the new code to reproduce 
exactly the results of Ref. 8 using the same physical 
model and chemical reaction rates. Specifically, the 
DSMC code employs constant rotational and vibra- 
tional collision numbers of 5 and 50, and the chemical 
reaction rates from Ref. 8 with the degree of vibration- 
dissociation coupling set to zero. The second subsec- 
tion investigates the effects of employing the improved 
physical models. For this purpose, the DSMC code 
employs variable rotational and vibrational collision 
numbers, vibrationally-favored dissociation, and the 
new chemical rate constants derived from the contin- 
uum expressions. In the third sub-section, comparison 
is made between the DSMC code using the new models 
and the continuum computations of Ref. 1. 

Comparison with Previous Results 
Using Old DSMC Models 

The new code is first assessed through comparison 
with DSMC computations presented by Moss ti all 8 
for the AFE at an altitude of 78 km (the flow condi- 
tions are given in Table 1). The free-stream tempera- 
ture and Knudsen number were 188 K and 1.2xl0“ 3 re- 
spectively. For compatability with the study reported 
m Ref. 8, constant rotational and vibrational collision 
numbers of 5 and 50 are employed. Also, the chemi- 
cal reaction rates from Ref. 8 are used, the degree of 
vibration-dissociation coupling is set to zero, and the 
shock standoff distance is specified as 0.110 m from 
the body. A total of 1,000 computational cells and 
over 100,000 simulated particles are employed in each 
of the calculations reported in the current work. 

The temperatures computed in the present study 
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(lines) are compared with data taken from Ref. 8 (sym- 
bols) in Fig. 1. It is clear that significant thermal 
nonequilibrium effects are present. In comparison to 
the profiles from Ref. 8, the present temperatures are 
slightly higher. This effect is at least partially due to 
the absence of radiative cooling in the current work. 

The mole fractions of the major species computed in 
the present study (lines) are compared with data taken 
from Ref. 8 (symbols) in Fig. 2. The present study is in 
very good agreement with the solutions from Ref. 8. In 
addition to the general form of the computed profiles 
shown in Figs. 1 and 2, the peak electron temperature 
of about 18,000 K and the maximum electron mole 
fraction of about 0.013 are in good correspondence to 
the results of Ref. 8. 

The radiative emission along the stagnation stream- 
line given by NEQAIR for the present DSMC flowfield 
solution is shown in Fig. 3. This profile compares quite 
well with the profile given in Ref. 8, although the peak 
value is somewhat larger in the present study. Integra- 
tion of the spontaneous emission along the streamline 
to the body gives a one-dimensional radiative flux of 
64.8 kW/m 2 /sr, which is about a factor of 1.6 greater 
than the result of Moss ti al. The difference may 
partly be accounted for in terms of radiative cooling 
which is omitted in the present calculations by estimat- 
ing radiative emission decoupled from the flowfield. 
Considering the significant differences in the radiation 
models employed in the two studies, the agreement is 
acceptable. 

The radiative heating at the stagnation point due to 
the radiative flux given above is found using a spherical 
cap model 1 . First, the infinite slab result is obtained 
for the total radiative heating flux at the stagnation 
point by multiplying the one-dimensional value by 2z 
for the optically thin gas part of the spectrum, and by 
x for the self-absorbing regions. The spherical cap re- 
sult is then computed as about 80% of the infinite slab 
value 1 . This procedure gives a total radiative heating 
value of about 340 kW/m 2 , which is a very high value. 
This value is even higher than the convective heating 
rate reported in Ref. 8 as 248 kW/m 2 . 

The contributions to the total radiative heating 
made by the various spectral regions which are ana- 
lyzed are listed in Table 3 under Old Model. The rea- 
sonable agreement found between the present DSMC 


results and those obtained in Ref. 8 for the radiative 
flux, indicates that the decoupled approach to radia- 
tion assumed in this study is an appropriate solution 
method for these flow conditions. In addition, analysis 
of the rate of loss of energy from the flowfield showed 
that only about 0.2% was due to radiative emission. 
This confirms that simulation of radiation coupled to 
the fluid mechanics is unnecessary for these flowfields. 

Comparison of D id and New DSMC 

Having established that the new DSMC code pro- 
vides solutions which are similar to previous studies, 
when using similar models, it is appropriate to assess 
the effect on the radiative emission by the introduc- 
tion of the improved physical models. Specifically, 
variable rotational and vibrational collision numbers, 
vibrationally-favored dissociation, and the new chemi- 
cal rate constants derived from the continuum expres- 
sions are implemented. The translational and vibra- 
tional temperatures computed with the DSMC code 
using the old and new models are shown in Fig. 4 for 
the same conditions considered in the previous subsec- 
tion. The relaxation zone is much larger with the im- 
proved models. Note that the vibrational mode using 
the new model does not equilibrate with the transla- 
tional mode until the body surface is reached. This is 
due mainly to the use of the variable vibrational col- 
lision number. This quantity is a strong function of 
temperature, and only reaches a maximum of about 
50 at the peak translational temperature. For most 
of the stagnation streamline, the vibrational collision 
number is greater than 100 which reduces significantly 
the rate of equilibration of the vibrational mode. 

The species which radiate most energy in the ultra- 
violet region of the spectrum (0.11-0.18^ m) are atomic 
nitrogen and oxygen. The variation in the mole frac- 
tions for these species is shown in Fig. 5. With the 
new, improved models, the rise in N and O due to 
dissociation is retarded significantly due to the use of 
the vibrationally-favored dissociation models and the 
lower vibrational temperatures. The strongest radi- 
ator above 0.2 /i m is N 2 + and its variation in mole 
fraction computed with the old and new DSMC mod- 
els is compared in Fig. 6. The rise in N 2 + is faster 
with the old models although behind the shock the 
two solutions show general agreement. 

The profiles of radiative emission are compared in 


AIAA 92-2971 5 


Fig. 7. The peak emission obtained with the new 
models is significantly lower than that computed with 
the previous models due to the lower electronic tem- 
perature, which is set equal to the vibrational value 
in NEQAIR. The computed three-dimensional inte- 
grated radiative heating at the stagnation point is only 
about 23.4 kW/m 2 for the new model which is a fac- 
tor of about 15 smaller than that predicted by the 
old model. The contributions to the total radiative 
heating made by the various spectral regions which 
are analyzed using NEQAIR are listed in Table 3 un- 
der New Model. The very large difference obtained in 
the radiative heating estimates using the old and new 
thermochemistry models demonstrates the sensitivity 
of the calculations to the physical models. 

Compa rison with Continuum Results 

The DSMC code using the new models is applied to 
the flow conditions examined by Whiting and Park 1 , 
which are given in Table 1. These are slightly different 
from those investigated in Ref. 8. The most important 
difference involves the shock standoff distance which is 
significantly larger in Ref. 1 (0.188 m) than in Ref. 8 
(0.110 m). The freestream temperature is 188 K and 
the Knudsen number is 1.7xl0~ 3 . The DSMC compu- 
tations again employed the energy-dependent proba- 
bilities of rotational and vibrational energy exchange, 
and the vibrationaily-favored dissociation model dis- 
cussed by Boyd. 4 The vibrational and chemical rate 
constants are made consistent with those employed in 
the continuum simulation as discussed earlier. 

The continuum results taken from Ref. 1 were com- 
puted using SPRAP which was developed by Park. 3 
This code computes the viscous, one-dimensional, con- 
tinuum flow behind an infinitely thin normal shock 
wave using the approach described in Ref. 3. The so- 
lution for the one-dimensional, uniform area duct may 
then be transformed to represent the flow along the 
stagnation streamline of a spherical body. The initial 
flow conditions immediately behind the shock are de- 
termined from the Rankine-Hugoniot shock-jump rela- 
tionships, assuming that the vibrational temperature 
remains at the free-stream value. A viscous shock layer 
treatment is applied for computation of the boundary 
layer which forms next to the cold wall of the vehicle. 

It is assumed that the flow in the boundary layer is 
chemically frozen with a ratio of specific heats equal 


to 11/9. Further, the wall is assumed to be noncat- 
alytic for the dissociation reactions and fully catalytic 
for the ionization reactions. In this approach, a two- 
temperature dissociation model is used in generating 
the reacting flow behind the normal shock wave. This 
model equates the molecular rotational temperature to 
the kinetic temperature of the atoms and molecules, 
and also equates the electron kinetic and electronic 
temperatures to the molecular vibrational tempera- 
ture. In the two-temperature model, all molecules 
have the same vibrational temperature, the degree of 
ionization is small, and the chemical reactions occur 
in the ground states of the atoms and molecules. 

The computed temperature profiles are compared 
in Fig. 8. The two- temperature continuum approach 
provides a translation-rotation value, and a vibration- 
electron-electronic value. The particle DSMC ap- 
proach provides separate values for the translational, 
rotational, and vibrational energy modes and equates 
the electronic temperature to the vibrational temper- 
ature to calculate the radiation. 

The temperature profiles exhibit many differences. 
The DSMC solution shows considerable shock thick- 
ness and structure which is omitted in the continuum 
calculation. In the DSMC computation, there is a sig- 
nificant region of the shock where the translational 
and rotational modes are not equilibrated, thereby 
casting doubt on the validity of the continuum two- 
temperature approach at least at these low densities. 

The rise in vibrational temperature predicted by 
DSMC is slower than that computed in the contin- 
uum simulation. Further, the DSMC method pre- 
dicts a higher degree of nonequilibrium between the 
vibrational mode and the translational and rotational 
modes all along the stagnation streamline. The rea- 
sons for these differences in the two solutions are not 
clear at this point. The factors involved include dif- 
ferences in simulating viscosity, thermochemistry, and 
transforming a one-dimensional calculation into a stag- 
nation streamline flow. These factors need to be in- 
vestigated more thoroughly. 

The mole fractions of atomic nitrogen and oxygen 
computed with the continuum method and the DSMC 
code for these conditions are compared in Fig. 9. 
Generally, the agreement is quite good except within 
the shock wave itself which is omitted in the con- 
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tinuum code. In the continuum technique, it is as- 
sumed that no chemical reactions occur upstream of 
the shock standoff location. The DSMC computation 
indicates that significant chemical activity does take 
place through the shock. This is due to the finite 
thickness of the shock-wave and to diffusion effects. 
Thus, the rise in atomic mole fractions computed by 
DSMC preceeds the continuum results by a significant 
distance. The profiles computed by the two methods 
then follow similar paths until the last few centimeters 
next to the body, where the continuum code assumes 
the flow to be chemically frozen. The variation in mole 
fraction of N 2 + computed with the two techniques is 
shown in Fig. 10. The rise in N 2 + predicted by DSMC 
occurs in the thick shock-front and therefore leads the 
continuum solution. Behind the shock, the DSMC re- 
sults provide a higher concentration of than the 
continuum data. 

The spontaneous radiative emission profiles ob- 
tained from the continuum and DSMC solutions are 
compared in Fig. 11. The peak emission obtained 
with the new DSMC models is significantly lower than 
that computed with the continuum method, due to the 
smaller electronic (vibrational) temperature. 

The integrated radiative heating at the stagnation 
point computed with the DSMC calculation is about 
47 kW/m 2 while the continuum solution gives a value 
of 75.5 kW/m 2 . The contributions to the total radia- 
tive heating made by the various spectral regions are 
listed in Table 4. 

Despite the large differences observed in the flow- 
field solutions, it is interesting to note that the con- 
tinuum and DSMC total radiative heating estimates 
lie within a factor of 2 of each other. It is also worth 
noting that by increasing the shock standoff distance 
from 0.110 m (from Ref. 8) to 0.188 m (from Ref. 1) 
the total radiative heating predicted by the new DSMC 
thermochemistry models is doubled, indicating nearly 
linear scaling between these two quantities. 

For the sake of completeness, the flow conditions 
investigated in Ref. 1 were also computed with the 
old DSMC thermochemistry models. The difference 
in flow quantities computed with the old and new 
models was similar to those shown in Figs. 4 through 
6. The flow quantities were again interpretted in 
terms of total radiative heating at the stagnation point 


using NEQAIR. The solution obtained with the old 
DSMC thermochemistry models gave a radiative heat- 
ing which was ten times higher than that predicted 
with the new models. 

Concluding Remarks 

A new approach has been evaluated for predict- 
ing the radiative heating of a blunt-body entering the 
Earth’s atmosphere. In this approach, the fluid me- 
chanics of the flow along the stagnation streamline was 
predicted using a particle method (the direct simula- 
tion Monte Carlo method, DSMC), which is decoupled 
from the radiative emission. 

Comparison was made with a previous DSMC cal- 
culation which employed outdated thermochemistry 
models together with a simplified, coupled radiation 
model. For the purpose of this comparison, the present 
DSMC computation also employed the old DSMC 
thermochemistry models. These two very different ap- 
proaches gave agreement to within a factor of 2 for the 
total radiative heating at the stagnation point. This 
is viewed as acceptable considering the differences be- 
tween the radiation models employed. The decoupled 
approach gave the higher value which was attributed 
in part to the absence of radiative cooling. 

For the same flow conditions, the use of new DSMC 
thermochemistry models led to a decrease by a factor 
of 15 in the total radiative heating at the stagnation 
point. This drastic reduction in heating illustrates the 
sensitivity of the computed data to the physical models 
employed in the simulations. It is proposed that the 
new thermochemistry models, which have been suc- 
cessfully evaluated in previous studies, should provide 
more realistic simulation results. This decrease in ra- 
diative heating indicates that it is not necessary to cou- 
ple radiation to the fluid mechanics under these flow 
conditions. To attain greater confidence in the numer- 
ical simulations, it is clear that experimental data is 
required for the validation of these physical phenom- 
ena. 

A further comparison of the new DSMC chemistry 
models with a continuum calculation gave agreement 
within a factor of 2 for the total radiative heating, 
despite considerable differences in the flowfield solu- 
tions. The main conclusion of this study is that there 
remains a large degree of uncertainty in the applica- 
tion of state-of-the-art numerical methods for the pre- 
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diction of radiative heating to hypersonic vehicles fly- 
ing in low-density regions of the Earth’s atmosphere. 
Considering the importance of such heating to many 
future aerospace projects, there is a need for continued 
numerical and experimental investigation in this area. 
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Table 1. Freestream conditions. 


Source Altitude Uoo (m/s) 

(kg/m 3 ) 

Ref. 8 78 km 

9110 

2.75x10-* 

Ref. 1 80 km 

9756 

2.05x10-* 

Table 2. Spectral regions analyzed with NEQAIR. 

Region (/im) 

Description 

Absorption 

0.1130-0.1140 

N atomic lines 

Yes 

0.1160-0.1170 

N atomic lines 

Yes 

0.1170-0.1180 

N atomic lines 

Yes 

0.1195-0.1205 

N atomic lines 

Yes 

0.1240-0.1250 

N atomic lines 

Yes 

0.1294-0.1314 N and 0 atomic lines Yes 

0.1315-0.1325 

N atomic lines 

Yes 

0.1323-0.1333 

N atomic lines 

Yes 

0.1405-0.1415 

N atomic lines 

Yes 

0.1490-0.1500 

N atomic lines 

Yes 

0.1740-0.1750 

N atomic lines 

Yes 

0.2000-2.0000 

Molecular lines 

No 

Table 3. Total radiative heating (in kW/m 2 ) using 
DSMC for the flow conditions from Ref. 8. 

Region (/xm) 

Old Model 

New Model 

0.1130-0.1140 

3.11 

0.08 

0.1160-0.1170 

4.19 

0.16 

0.1170-0.1180 

2.30 

0.08 

0.1195-0.1205 

6.07 

0.21 

0.1240-0.1250 

3.65 

0.14 

0.1294-0.1314 

8.23 

0.24 

0.1315-0.1325 

3.16 

0.23 

0.1323-0.1333 

2.37 

0.07 

0.1405-0.1415 

3.13 

0.28 

0.1490-0.1500 

11.94 

0.72 

0.1740-0.1750 

19.30 

2.21 

Sub- Total 

67.5 

4.42 

0.2000-2.0000 

272.2 

19.0 

Total 

339.7 

23.4 


Table 4. Continuum and DSMC contributions (in 
kW/m 2 ) of the various spectral regions to the total 
radiative heating at the stagnation point for the flow 


conditions from Ref. 1 at an altitude of 80 km. 


Region (/im) 

Continuum 

New Model 

0.1130-0.1140 

0.44 

0.20 

0.1160-0.1170 

0.76 

0.50 

0.1170-0.1180 

0.32 

0.27 

0.1195-0.1205 

0.84 

0.46 

0.1240-0.1250 

0.92 

0.40 

0.1294-0.1314 

0.76 

1.13 

0.1315-0.1325 

0.76 

0.61 

0.1323-0.1333 

0.52 

0.19 

0.1405-0.1415 

0.84 

0.70 

0.1490-0.1500 

2.84 

1.65 

0.1740-0.1750 

5.36 

4.73 

Sub- Total 

14.4 

10.8 

0.2000-2.0000 

61.1 

36.2 

Total 

75.5 

47.0 
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Fig. 1. Temperature profiles along the stagnation 
streamline: comparison of present calculator and 
Ref. 8 using the old DSMC thermochemistry models. 



Fig. 2. Profiles of atomic species mole fractions: com- 
parison of present calculator and Ref. 8 using the old 
DSMC thermochemistry models. 



Fig. 3. Profiles of radiative emission along the stagna- 
tion streamline: comparison of present calculator and 
Ref. 8 using the old DSMC thermochemistry models. 



o 5 10 15 20 


Distance from body (cm) 

Fig. 4. Temperature profiles along the stagnation 
streamline: comparison of the old and new DSMC 
thermochemistry models. 



Distance from body (cm) 

Fig. 5. Profiles of atomic species mole fractions along 
the stagnation streamline: comparison of the old and 
new DSMC thermochemistry models. 
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Fig. 6. Profiles of Ns + mole fraction along the stagna- 
tion streamline: comparison of the old and new DSMC 
thermochemistry models. 



Fig. 7. Profiles of radiative emission along the stagna- 
tion streamline: comparison of the old and new DSMC 
thermochemistry models. 



Fig. 8. Temperature profiles along the stagnation 
streamline: comparison of continuum and DSMC cal- 
culations using the new thermochemistry models. 



Fig. 9. Profiles of atomic species mole fractions: com- 
parison of continuum and DSMC calculations using 
the new thermochemistry models. 



Fig. 10. Profiles of N 2 + mole fractions: comparison 
of continuum and DSMC calculations using the new 
thermochemistry models. 



Distance from body (cm) 

Fig. 11. Profiles of radiative emission: comparison 
of continuum and DSMC calculations using the new 
thermochemistry models. 
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SIMULATED PITCH, YAW, AND ROLL TORQUES 
ON THE MAGELLAN SPACECRAFT DURING AEROBRAKING 

Brian L. Haas* 

Eloret Institute , 3788 Fabian Way , Palo Alto , California 94303 
Durwin Schmittf 

Martin Marietta Astronautics Corp,, Denver, CO 

Aerodynamic loads upon the Magellan spacecraft during proposed acrobraking through the atmo- 
sphere of Venus are computed at off-design attitudes with a direct Monte Carlo particle simulation 
method \ This method is not restricted to the assumption of collisionless flow normally employed to 
assess spacecraft aerodynamics . Simulated rarefied flows at the nominal altitude of 140 km and entry 
speed of 8.6 km/s were compared to simulated and analytic free molecular results. The nominal space- 
craft orientation is defined with the flow directed along the central body axis and perpendicular to the 
solar panels. Aerodynamic moments and lift for rarefied entry at several angles of pitch and yaw were 
at least 10% greater than free molecular results, while drag was about 20% lower. At each orientation , 
the resulting aerodynamic center of pressure on the body axis was located behind the vehicle center of 
gravity ; providing restoring torques and promoting spacecraft stability. Roll torques about the body axis 
for entry at nominal orientation , but with the solar panels canted at various angles to the freestream flow ; 
were computed and compared to free molecular results. This configuration has been suggested as an inno- 
vative experiment to assess gas-surface interaction during acrobraking. Periodic free-molecule boundary 
conditions and a coarse computational grid and body resolution served to minimize the simulation size 
and cost while maintaining solution accuracy. 


NOMENCLATURE 

n number density 

T temperature 

u flow velocity 

a intermolecular potential exponent 

e surface radiative emissivity 

A gas mean free path 

p gas density 

Subscripts: 

s stagnation value 

oo freestream value 


INTRODUCTION 

The Magellan spacecraft has been mapping the surface 
of Venus from a highly elliptic orbit (eccentricity, e = 0.39) 
since September of 1990. Mission planners at NASA’s Jet 
Propulsion Laboratory (JPL) would like to circularize the 
orbit to improve mapping but cannot perform the manuever 
through thruster activity alone due to limited remaining pro- 
pellant. The orbital maneuver may be achieved through a 
series of gentle passes through the Venus atmosphere (1600 
passesatu = — 1 m/s each), 1 Besides circularizing the orbit, 
these maneuvers would provide considerable data related to 
atmospheric entry of satellite spacecraft. 


* Research Scientist. Member, AIAA. 

Mailing Address: NASA Ames Research Center 
M/S 230-2, Moffett Field, California 94035-1000 
f Retired. Member, AIAA. 

This paper ii declared a work of the U.S. Government and is not 
subject u> copyright protection in the United S utes. 


The spacecraft configuration is depicted in Fig. 1 in its 
nominal entry orientation with flow directed (at orbit pe- 
riapsis) along the central body axis. The solar panels in 
this configuration are normal to the flow, but indeed may 
be canted at any angle to form an effective “windmill” out 
of the spacecraft during entry. Restoring roll torques on the 
spacecraft may be measured along with surface heating and 
orbit altitudes to deduce the flowfield density and the sur- 
face accomodation coefficients in the normal and tangential 
directions. Axes for pitch and yaw are indicated in the front 
view. 

Aerodynamic heating on the delicate solar panels dur- 
ing each aeropass was one mission constraint. Direct parti- 
cle simulations of entry at several altitudes from 125 km to 
140 km in the nominal orientation 2 verified that this heat- 
ing was within the specified tolerance at altitudes exceeding 
135 km. At the nominal altitude of 140 km, the heating and 
drag results corresponding to the rarefied flow were very 
close to those corresponding to free molecular (collision- 
less) flow at that altitude. 

Mission planners are now concerned with spacecraft sta- 
bility during acrobraking, particularly if the entry orienta- 
tion has a high angle of pitch or yaw. Simplified analyses 
rely upon the assumption of free molecular flow which may 
not be valid at nominal flight conditions. The objective of 
the present study is to calculate body torques on the Magel- 
lan spacecraft for several entry orientations at the nominal 
entry altitude and periapsis velocity (u^ = 8.6 km/s), em- 
ploying a particle simulation method which is not restricted 
to the free molecule flow assumption. The results of the 
rarefied simulation may then be compared to free molecu- 
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lar results to assess the significance of molecular collisions 
in the flow. Roll torques associated with the windmill entry 
configurations were also computed at nominal conditions 
and compared to free molecule results. The roll results are 
not presented In this abstract but will appear in the confer- 
ence paper. 

PARTICLE SIMULATION METHOD 

Direct simulation Monte Carlo (DSMC) particle meth- 
ods model the motion and interaction of thousands or mil- 
lions of computational particles to simulate gas dynamics. 3 
Specifics of the DSMC code employed in these studies are 
provided in Ref. 2. Given a particular position, velocity, 
and internal energy status, each particle in the flowfield 
travels unobstructed along the linear trajectory of its veloc- 
ity vector over the duration of a single time step. At this 
time, neighboring particles are identified throughout the 
flowfield and paired-off as potential collision candidates. 
The flowfield is divided into a network of cells to facil- 
itate identification of neighboring particles and to define 
the finest resolution for sampling macroscopic (low prop- 
erties. Employing probabilities as functions of individual 
collision parameters such as collision cross-section and rel- 
ative translational speed 4 , the subset of all candidate pairs 
which collide during the timestep are identified. In simulat- 
ing free molecule flow, the probability of collision is arti- 
ficially set to zero, representing an infinite molecular mean 
free path. 

The entire simulated flowfield is initialized with free 
stream conditions. The particle simulation then runs 
through a transient phase as the solution develops and flow- 
field structures form. Upon reaching steady-state, the sim- 
ulation collects statistical samples for measuring properties 
of the flowfield and body surfaces. 

The code employed in the present study was developed 
by McDonald 5 for efficient implementation on vector su- 
percomputers. This code simulates non-reactive 3D (low 
of general gas mixtures about arbitrary geometries. Molec- 
ular collisions pertain to the Variable Hard Sphere model 
of Bird 6,4 with an inverse intermolecular potential expo- 
nent of a = 5. The internal energy modes are modeled 
with three fully-excited degrees of freedom to account for 
molecular rotation and vibration. Internal energy excitation 
is performed with the mechanics of Borgnakke and Larsen 7 
employing a fixed probability of relaxation of 1/5. 

Body Geometry and Grid 

The geometry of the Magellan spacecraft is shown in 
Fig. 1 and compared to the simulated geometry. Due to 
the large mean free path in the flow about the spacecraft, 
small features on the vehicle such as the altimeter antenna 
(ALTA), the medium and low gain antennas (MGAJ.GA), 
and the rocket engine modules (REM) have negligible im- 
pact on the flowfield as a whole and may be excluded to 
simplify the simulation geometry, leading to two planes of 
symmetry on the vehicle. 


To represent a surface in the cubic cartesian grid network 
of the simulation, it is necessary to approximate the surface 
as a composite of planar facets. Each facet has a normal de- 
fined from the intersections of the surface with the edges of 
the cell. This faceted description of the body is appropriate 
given that body radii are large in comparison to the cell size, 
and that the intersection of different surfaces occurs at cell 
boundaries. Since the solar panels are a dominant compo- 
nent of the structure, the simulation cell network was scaled 
such that the square simulated panels (measuring 8 cells x 8 
cells) have the same frontal area as the actual solar panels. 
This leads to the scaling factor, 1 cell = 31.44 cm. Such a 
coarse grid and body resolution was used in this study to 
reduce the computational expense associated with running 
several simulations. 

Grid resolution greatly impacts the accuracy of solutions 
for vehicle aerodynamics and heating in rarefied flows. An 
established criterion for sufficient grid resolution is that the 
local mean free path must exceed the cell dimension. For 
cold-wall blunt-body rarefied (lows, flowfield density will 
likely be quite large near the body surface, leading to small 
stagnation mean free paths. Simulated at nominal flight 
conditions, the (low density n and temperature T along the 
stagnation streamline at the center of the solar panel are 
plotted in Figs. 3 and 4. Equilibrium kinetic theory pro- 
vides the following estimate of the local mean free path 3 at 
the body surface. 


At nominal conditions 9 the freestream mean free path 
' s = 23.4 m (74.43 cells). The stagnation properties 
lead to a stagnation mean free path around A s % 4.2 cells, 
exceeding the minimum accuracy criterion. 

Body surfaces are modeled as if in radiative equilibrium 
with deep space at a temperature of 4 K and emissivity 
of e = 0.82, and employ a thermal accommodation coeffi- 
cient of unity, typical of rough cool surfaces facing into the 
flow. This model is simple to implement in the simulation, 
it does not require a prescribed estimate of surface temper- 
ature, and it allows each surface facet to reach its own tem- 
perature independent of neighboring facets. It is assumed 
that radiation from the flowfield or from other body surfaces 
would contribute negligibly to the net heating of a given 
surface facet. 

Flow Boundary Conditions 

In simulating highly-rarefied flows, the computational 
flow domain must extend far enough upstream of the body 
to provide ample opportunity for freestream molecules to 
interact with those molecules that have reflected from the 
body and are diffusing upstream. Insufficient upstream do- 
main size leads to overprediction of aerodynamic heating 
and forces. 

Taking advantage of body symmetry, the simulated flow- 
field configuration is that of a wind tunnel depicted for pitch 
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simulations in Fig. 2 (100 L x 120 Hx 12 Win cells) with a 
specularly-reflecting symmetry plane and a free-molecular 
outer plane. Particles which strike the outer plane after re- 
flecting off of a body surface are removed from the flow 
field. All other particles striking the outer plane reflect 
specularly back into the flow. This effectively reduces the 
influence that the outer wall has upon the flow near the ve- 
hicle in highly-rarefied or free molecular flows. 

Particles enter the flow domain from the left plane at 
the specified pitch angle, and exit normally from the right 
plane. The top and bottom planes represent free-molecular 
periodic boundary conditions. Again, particles which had 
reflected from the body before striking these planes are re- 
moved from the simulation. All other particles striking 
one of these planes simply re-enter the flow domain from 
the opposite plane, maintaining its original velocity vector. 
This boundary condition still provides freestrcam flow ap- 
proaching the body geometry while reducing the required 
height of the wind tunnel, therefore minimizing the size and 
cost of the simulation. 

SIMULATION RESULTS 

Employing a freestream flow speed of = 8.6 km/s 
and the day-side atmospheric data of Keating 9 at 140 km 
altitude, simulations were performed for several entry ori- 
entations of pitch and yaw. Results from the particle sim- 
ulation for rarefied flows, employing the finite molecular 
mean free path of A^ = 74.43 cells, were compared to free 
molecular (collisionless) simulation results. These were, in 
turn, compared to analytic free molecular results generated 
with the FREEMAC code 10 which employs ray-tracing al- 
gorithms to assess aerodynamic coefficients of spacecraft. 

Limited computational resources and the large number 
of cases investigated in the present study restricted the num- 
ber of particles which could be employed in the simulauon 
to just four particles per cell in the freestream. However, 
density gradients in the flowfield significantly increased the 
particle densities near body surfaces such that 400,000 par- 
ticles existed in the flowfield at steady-state. Employing 
roughly 4,000 transient steps and 6,000 steady-state sam- 
pling steps in the simulation, the number of statistical sam- 
ples was sufficiently large to yield accurate and meaning so- 
lutions. Run-times averaged 1. 0-2.0 CPU hours per case 
on the Cray-YMP supercomputer. The simulation could 
have employed ten times as many particles but would have 
required a ten-fold increase in computational time per case 
which was not warranted in the present study. 

The simulation computes the net force and heat flux 
upon each surface facet of the vehicle. Torques are com- 
puted by the moments of body forces about a reference 
point defined by the intersection of the central body axis 
and the solar panel axis. This reference point is very close 
to the spacecraft center of mass. Moment coefficients are 
defined using the frontal area of the simulated geometry 
(Aj = 236.8 square cells), a reference length given by 
the HGA diameter ( D = 11.77 cells), the freestream den- 


sit y (Pco = 6. 13 x 10~ 9 kg/m 3 ), and the freestream velocity. 

The moment coefficients at each pitch attitude are given 
in Fig. 5. Note that the free molecular results from the par- 
ticle simulation agree very well with FREEMAC results. 
Moment coefficients in the rarefied flow simulation, how- 
ever, differ from free molecular results by as much as 40%. 
For all cases, the moment was directed so as to reduce the 
pitch of the vehicle. This was also observed in the yaw atti- 
tude cases plotted in Fig. 6. Again there is good agreement 
between the free molecule simulation and FREEMAC. Rar- 
efied yaw results were not available at the time of writing 
this abstract, but will appear in the conference paper. 

Aerodynamic lift and drag are plotted in Fig. 7 for 
the pitch cases. The limited data points available from 
FREEMAC agree very well with the particle simulation re- 
sults for drag in free molecular flow, but flowfield collisions 
in the rarefied flow led to lower drag. Pitching the vehicle 
led to a downward force (negative lift) which was nearly 
twice as great in the rarefied flows than in the collisionless 
flows. 

Finally, the location Xcp of the center of pressure be- 
hind the moment reference point on the central body axis 
is plotted in Fig. 8 for both pitch and yaw cases. At all 
attitudes, x^ lies behind the center of gravity, promoting 
spacecraft static stability. Again, molecular collisions in 
the rarefied flows do have a noticeable effect upon x^. 

The con Terence paper will include a comparison between 
FREEMAC results for the coarse geometry as used in the 
particle simulation, and a very fine geometry which closely 
resembles the actual spacecraft configuration. These results 
can be generated very quickly, but were not a vailable at the 
time of writing this abstract. 

CONCLUDING REMARKS 

The aerodynamics of the Magellan spacecraft during 
proposed entry into the atmosphere of Venus may be com- 
puted for free molecular flows with the FREEMAC code. 
However, at the nominal altitude of 140 km, the free molec- 
ular flow assumption may no longer be valid. This study 
computed ihe body torques on the spacecraft when enter- 
ing with various-off-design attitudes, using a direct particle 
simulation method for rarefied flow conditions. Repeating 
these simulations for collisionless flow provided a direct 
means of assessing the significance of molecular interac- 
tions in the flow, as well as providing a means to validate the 
technique through comparison to FREEMAC results. Use 
of a coarse body resolution permitted repeated inexpensive 
calculations without sacrificing solution accuracy. Results 
indicate that body torques will act to restore the vehicle to 
its nominal zero-pitch, zero-yaw attitude. The quantitative 
results of this study may facilitate assessment of spacecraft 
stability. Most importantly, this study demonstrates that the 
assumption of free molecular flow underpredicts the result- 
ing body torques and lift, while overpredicting drag. 
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Fig. 1 Magellan spacecraft configuration versus simulated geometry shown at 
nominal attitude (zero pitch and yaw). 
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Fig. 6 Coefficient of restoring yaw moment for entry at various yaw angles. ^'8* 8 Location of the center of pressure along the central body axis behind 

Comparison of free molecular particle simulations results to FREEMAC re- reference point (approximate vehicle center of mass) for entry at various 

sl ,l ls pitch and yaw attitudes. Comparison of rarefied and free molecular particle 

simulations results for pilch. 
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Fig, 5 Results of simulation employing dual-node model (Model D): (a) transient heat- 7 Effects of heat transfer parameters c f r and k, 0 upon surface temperatures for 

lug due to convection Q, conduction Vk> and radiation q r ; and (h) transient tempera- ocrtibraklng using Model 1) at perigee conditions for 0=5.6°. Define 7, as the limit 
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FLOW RESOLUTION AND DOMAIN OF INFLUENCE 
IN RAREFIED HYPERSONIC BLUNT-BODY FLOWS 

Brian L. Haas* 

Eloret Institute, 3788 Fabian Way, Palo Alto, California 94303 


A direct particle simulation method is used in this parametric study to assess the influence of upstream 
domain size and grid resolution upon how properties and body aerodynamics m rarefied flo ws over cold 
blunt bodies Insufficient grid resolution leads to over prediction of aerodynamic heating and drag, and 
JJSSSta of aerodynamic force perpendicular to the freestream flow direction. Solution accuracy 
within one percent error is attained when the grid cell size is no larger than roughly two 
free paths near the body. Molecular diffusion from the body surface results in greater ar-fie d domain 
influence as flow rarefaction increases. As a result, insufficient upstream domain size leads to overpre- 
diction of heating and drag along with underprediction of perpendicular force. However, since the use o 

size are of less consequence as freestream rarefaction increases. Defining x Tm „ as the position ahead 
oTthTb^y where peak temperature occurs, this study shows that one percent error in solution accuracy 
is attained when the upstream computational domain size exceeds 3x Tm% .. Simulation of a hard-sphe e 
giTmore sensitive to grid resolution while simulation of a Maxwell gas is more sensitive to upstream 

domain size . 


NOMENCLATURE 

£} cylinder diameter 

F aerodynamic force 

Kn Knudsen Number 

l ! upstream domain size 

M Mach number 

n number density 

q aerodynamic total heat flux 

T temperature 

x location relative to stagnation point 

a intermolecular potential exponent 

A gas mean free path 

Subscripts: 

s stagnation value 

w wall value 

x in flow direction, drag 

y perpendicular to flow direction 

co freestream value 

INTRODUCTION 

Flow field characteristics about blunt bodies during at- 
mospheric entry lead to considerable challenges for com- 
putational simulation. Highly rarefied flows, where the gas 
density n is low and the mean free path A between molecu- 
lar collisions is high, are better suited to particle simula- 
tion methods, such as the direct simulation Monte Carlo 
(DSMC) technique pioneered by Bird, 1 rather than contin- 
uum techniques based upon the Navier-Stokes equations. 


* Research Scientist. Member, AIAA. 

Mailing Address: NASA Ames Research Center 
M/S 230-2, Moffett Field, California 94035 

This paper is declared a work of ihe U.S. Government and is not 
subject to copyright protection in the United States. 


DSMC methods employ many model particles whose mo- 
tion and interaction simulate gas dynamics directly. The 
simulated flowfield is divided into a network of small cells 
to facilitate collision modeling and statistical sampling. 
The computational burden of DSMC methods, however, 
grows proportionally with local gas density and the size of 
the computational domain. In typical entry flows, the body 
surface temperature is rather low compared to the stagna- 
tion temperature, leading to a steep density gradient near 
the body surface. Accurate simulation of this flow requires 
sufficient grid resolution near the body. Alternatively, the 
extent of freestream rarefaction results in a leading shock 
layer which is fully merged with the boundary layer of 
the body, yet extends far upstream. Accurate simulation 
therefore requires a large upstream computational flow do- 
main. Together, these requirements rapidly drive the com- 
putational cost of the particle simulation method upward. 
The researcher must therefore understand what simulation 
costs are necessary to obtain a solution of sufficient accu- 
racy. The objective of this parameteric study is to assess 
quantitatively the sensitivity of aerodynamic loads and gas 
properties to grid resolution and simulation domain. 

The present study employs an efficient particle simula- 
tion technique 2,3 to investigate the net heat flux q and aero- 
dynamic forces F on two-dimensional circular cylinders of 
diameter D, along with gas temperature T and number den- 
sity n along the stagnation streamline, for hypersonic flows 
at Mach number M ^ = 20 over the range of Knudsen 
numbers given by Kn= A/D = {0.1,0.3,1.0,3.0,10.0}. 
The gas is modeled by hard-sphere particles with two fixed 
internal degrees of freedom. Interaction of the gas with 
the cylinder surface at wall temperature T w /T s = 0.05 
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is modeled as fully diffuse and thermally accommodated, 
where T 9 is the stagnation (total) temperature of the flow. 
All simulations employed between 10 and 20 particles per 
cell as the input freestream number density, n^. For each 
set of conditions, the simulation resolution and domain 
were varied considerably and their effects upon solution ac- 
curacy were determined. Although not presented in this ab- 
stract, results for the conference paper will include studies 
at = {5, 10} and T w /T c = {0.1, 0.5, 1.0}. 

In the DSMC technique, any two particles in the flow 
may collide only if they both reside in the same compu- 
tational cell. Likewise, any two particles within a given 
cell can collide regardless of their respective positions in 
the cell. If the cell dimension near the cold body surface is 
too large, then energetic particles at the far edge of the cell 
are able to readily transmit momentum and energy to par- 
ticles immediately adjacent to the surface. The latter par- 
ticles may, in turn, transmit that energy and momentum to 
the surface. This leads to over-prediction of both the surface 
heat flux and the aerodynamic force on the body than would 
occur in the real gas. This error is minimized by reducing 
the cell dimension relative to the local mean free path of 
molecules near the surface. From kinetic theory, the local 
mean free path A is related to the local number density n 
and temperature T as follows, 



fir) (£ 


2/a 


( 1 ) 


where subscript oo denotes values at freestream conditions, 
and a is the exponent of the assumed inverse-power inter- 
molecular potential which may vary between the limits of 
the Maxwell molecule (a = 4) and the hard sphere (q = co). 

As a consequence of Eq. (1), regions of high density lead 
to short mean free paths, requiring finer cell resolution to 
capture flow gradients accurately. Density profiles along 
the stagnation streamline for the cylinder flows above are 
plotted in Fig. 1. Note that the density at the cylinder sur- 
face (at x/D = 0) is significantly higher than the freestream 
value. This behavior is more pronounced, and the gradi- 
ents are steeper, for flows at lower Kn^. This, combined 
with the fact that lower Kn^ means that the freestream A M 
is short, dictates that the cell resolution must be very fine 
compared to more rarefied flows. 

Temperature profiles along the stagnation streamline for 
these same flows are also plotted in Fig. 1. These curves 
clearly demonstrate how far upstream of the body its pres- 
ence affects the flow. This domain of influence increases 
with Kn^ as a result of upstream diffusion of particles 
which reflected from the body surface. 

Flow density and temperature at the stagnation point on 
the cylinder may be used in Eq. (1) to determine the stag- 
nation mean free path A,. As measured in units of cell- 
lengths, A, is dependent upon the cylinder diameters D em- 
ployed in the simulations and is plotted for each Kn^ flow 
in Fig. 2. 

The plot in Fig. 3 compares the local Knudsen number 
along the stagnation streamline for a Maxwell-molecule gas 


and a hard-sphere gas. The gradients for the hard-sphere are 
steeper and lead to lower Kn at the stagnation point on the 
body compared to the Maxwell-molecule. Alternatively, 
the domain of influence extends much further upstream for 
the Maxwell-molecule. As a consequence, accurate DSMC 
solutions for hard-sphere gases require finer cell resolution 
while solutions for Maxwell-molecule gases require larger 
upstream computational fields. 

For each of the flow conditions above, the grid resolution 
is defined by the cylinder diameter D while the flow domain 
size is defined by the distance L 1 of the computational field 
upstream of the cylinder. These length scales are depicted 
schematically in Fig. 4 along with r Tmxx which represents 
the position of the peak temperature along the stagnation 
streamline upstream of the cylinder. 

GRID RESOLUTION STUDY 

For rarefied flows in the lower Knudsen number range 
0.1<Kn<1.0, simulations were performed to assess the ef- 
fects of grid resolution upon solution accuracy. At each 
value of Kn^ , the flow was simulated using several grid- 
resolutions corresponding to cylinder diameters of D = 
{5, 10, 20, 30, 50, 100, 200} cells-lengths. For all cases, a 
sufficient upstream domain size was used (L 1 j D = 2). 

The density and temperature profiles along the stagna- 
tion streamline ahead of the body are plotted for Kn^ =0.1 
in Fig. 5. For the coarse grids, where D< 50 cells, the pro- 
files were smeared considerably. For D> 50, the profiles 
coalesced to a single form and are therefore assumed to rep- 
resent accurate solutions. This same general behavior was 
observed for the other Kn cases, although the resolution re- 
quired to capture the appropriate profiles did not need to be 
as fine as in the Kn^ =0.1 flow. Indeed, for Kn^ = 1.0, 
accurate profiles were obtained at D> 20. 

For each simulation, the aerodynamic heating q and 
forces F x and F y were computed. Here, F x represents drag 
force in the direction of the freestream flow. F y is the net 
force, perpendicular to the freestream flow, on the top half 
of the body only. Note that the net lift force, integrated over 
the top and bottom halves of any axisymetric body, would 
equate to zero. 

To assess the error corresponding to a given flow res- 
olution, the body aerodynamics computed with the finest 
resolution at each Kn value was assumed to represent the 
“correct” solution for that flow. Employing the stagnation 
density and temperature from that solution, the stagnation 
mean free path. A, is determined from Eq. (1) for each res- 
olution case and plotted in Fig. 2. Errors in body aerody- 
namics are plotted against A, for all cases in Fig. 6. At all 
resolutions, heating and drag were overpredicted while F y 
was underpredicted (the absolute value of its error is plot- 
ted in the figure). These errors correlate fairly well with A, . 
Note that errors in heating are worse than errors in the forces 
(which are roughly equal but opposite). These plots indi- 
cate that aerodynamic errors will be less that one percent 
for grids in which the stagnation mean free path is on the 
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order of a half-cell or greater, regardless of the freestream 
rarefaction Kn^ . This conclusion, summarized by 

A, > 0.5 cells -iResolution Criterion (2) 

represents an appropriate design-criterion for grid resolu- 
tion in cold-wall blunt-body simulations. 

FLOW DOMAIN STUDY 

For rarefied flows in the higher Knudsen number 
range (1.0<Kn<10.0), simulations were performed to as- 
sess the influence of the upstream flow field domain upon 
solution accuracy. In view of the results above, the cell res- 
olution used in each case here was sufficiently fine such that 
the stagnation mean free path was roughly one cell-length 
near the body. All cases employed the hard-sphere molec- 
ular model. The objective here was to employ different up- 
stream computational field sizes, defined by length L f in 
Fig. 4, to determine the minimum acceptable domain size 
required for a given penalty in solution accuracy. 

Row density and temperature along the stagnation 
streamline is plotted in Fig. 7 at Kn^ =1.0 employing dif- 
ferent domain sizes L l . With sufficiently large L 9 , all the 
curves coalesced to a single curve to represent the assumed 
correct profile as represented at L f / D = 4.0. With insuf- 
ficient upstream domain (i.e. L f is small), density is un- 
derpredicted while temperature is overpredicted. However, 
once the domain size that is used exceeds about twice the 
distance of the peak temperature location the tem- 

perature and density profiles followed the expected forms 
predicted when ample upstream domain was employed. 
This is observed in Fig. 7 at L f / D = 0.7. It is somewhat 
surprizing that, although the flow was not simulated fur- 
ther upstream, the profiles matched the correct form as if 
it were. Though not presented here, behavior for flows at 
other Kn^ are qualitatively similar to the results above. 

Quantitative assessment of the errors in aerodynamic 
loads on the cylinder for simulations employing different 
upstream domain sizes are presented in Fig. 8. The er- 
rors for each case are computed relative to the aerodyan- 
mic loads which were obtained using the largest domain L 1 
and presumably represent the correct values. Each plot em- 
ploys a length scale normalized by * TmM . As expected, 
the errors in predicted aerodynamics decrease as larger up- 
stream domain sizes are used. Aerodynamic heating q and 
drag force F x were overpredicted in all cases while the force 
perpendicular to the freestream flow direction, F y , was un- 
derpredicted. 

Note that each curve begins to level off as L f jx Ttatk% 
drops below roughly 0.3. In the limit of decreasing up- 
stream domain, the incoming flow is almost completely un- 
affected by particles reflecting from the cylinder surface 
such that the aerodynamic loads approach the values cor- 
responding to the free molecular flow limit For heating 
and drag, the limit is less erroneous for flows at greater 
freestream Knudsen numbers Kn^ since free molecular 
flow represents the limit of greatest possible rarefaction. 


The interesting paradox which has developed is that, while 
the domain of upstream influence increases with Kn^ , the 
computational penalty in terms of solution accuracy is less 
dramatic compared to simulations at lower Kn^ . Regard- 
less, an appropriate criterion for defining a suitable compu- 
tational flow domain for rarefied flow simulation is given 
by 

l f /x >3. Domain Criterion (3) 

i Tm»* 

CONCLUDING REMARKS 

When employing direct particle simulation methods, ac- 
curate simulation of highly rarefied flows about cold blunt 
bodies requires sufficient computational grid resolution and 
upstream flow domain size. The objectives of this para- 
metric study were to assess quantitatively what penalty is 
suffered in solution accuracy when minimizing computa- 
tional expense by using short domains and coarse grids. 
Qualitative results were similar for all freestream Knud- 
sen numbers. Solution errors on the order of one percent 
for aerodynamic heating and forces are attained with grids 
for which the molecular mean free path near the stagna- 
tion point exceeds a half-cell in length while the upstream 
domain size exceeds three times the distance to the peak 
temperature point along the stagnation streamline. Simu- 
lated hard-sphere gases are more sensitive to grid resolution 
than Maxwell gases, but require less extensive upstream do- 
mains. The results presented in this abstract required less 
than two weeks to generate and will be duplicated in the 
final paper to include solutions at other freestream Mach 
numbers and surface temperatures. 
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Fig. 1 Profiles along the stagnation streamline for flows of differing rarefaction at A/*, - 20 
over cold cylinders; (a) number density, n; (b) temperature, T. 



Fig. 2 Molecular mean free path at the stagnation point of the cylinder of diameter D (both 

measured in units of cell-lengths) for freestream Knudsen numbers, 0.1 < Kn M < 10 

at = 20 and T w /T, = 0.05. 
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Fig. 3 Local Knudsen number along the stagnation streamline upstream of a cold cylinder 
for hard-sphere and Maxwell gases at Kn^ =0.1. 



Fig 4 Definition of simulation length-scales and schematic representation of temperature 
and density profiles upstream of a cold cylinder of diameter D in rarefied flow. 
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Fig 5 Profiles along the stagnation streamline upstream of a cold cylinder for flows simu- 
lated at different grid resolutions for Kn^ =0.1; (a) number density, n; (b) temperature, T. 
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Fig. 6 Errors in aerodynamic loads upon a cold cylinder for flows simulated at different grid 
resolutions; (a) heating, 9 ; (b) drag force, F r ; (c) perpendicular force, F y . Absolute values 
of error fractions measured relative to loads computed at the finest resolutions. 
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Fig. 7 Profiles along the stagnation streamline upstream of a cold cylinder for flows simu- 
lated with different upstream domain sizes for Kn^ = 1.0; (a) number density, n; (b) tem- 
perature, T. 
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Fig. 8 Errors in aerodynamic loads upon a cold cylinder for flows simulated with differ- 
ent upstream domain sizes; (a) heating, q\ (b) drag force, F r ; (c) perpendicular force, F y . 
Absolute values of error fractions measured relative to loads computed with the largest do- 


mams. 






